diff --git a/Project.toml b/Project.toml index 0a22707..c56bceb 100644 --- a/Project.toml +++ b/Project.toml @@ -5,12 +5,14 @@ version = "1.3.0" [deps] CEnum = "fa961155-64e5-5f13-b03f-caf6b980ea82" CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298" +GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f" NetworkOptions = "ca575930-c2e3-43a9-ace4-1e988b2c1908" PROJ_jll = "58948b4f-47e0-5654-a9ad-f609743f8632" [compat] CEnum = "0.2, 0.3, 0.4" CoordinateTransformations = "0.6" +GeoFormatTypes = "0.4" PROJ_jll = "900.100" julia = "1.6" diff --git a/README.md b/README.md index 8544fbe..0e8c8fc 100644 --- a/README.md +++ b/README.md @@ -11,8 +11,10 @@ Quickstart, based on the [PROJ docs](https://proj.org/development/quickstart.htm using Proj # Proj.jl implements the CoordinateTransformations.jl API. -# A Proj.Transformation needs the source and target coordinate reference systems. +# A Proj.Transformation needs the source and target coordinate reference systems (CRS), +# or a single pipeline. trans = Proj.Transformation("EPSG:4326", "+proj=utm +zone=32 +datum=WGS84") +# The CRS can be a string or the CRS type, which also interfaces with GeoFormatTypes.jl. # Once created, you can call this object to transform points. # The result will be a tuple of Float64s, of length 2, 3 or 4 depending on the input length. diff --git a/gen/generator.jl b/gen/generator.jl index 37d936b..a6a95ae 100644 --- a/gen/generator.jl +++ b/gen/generator.jl @@ -24,6 +24,8 @@ function rewrite(ex::Expr) fargs′ = [:x, :y, kw(:z, 0.0), kw(:t, Inf)] elseif fname in (:proj_create_crs_to_crs, :proj_create_crs_to_crs_from_pj) fargs′[3] = kw(:area, :C_NULL) + elseif fname in (:proj_get_id_auth_name, :proj_get_id_code) + fargs′[2] = kw(:index, 0) elseif fname === :proj_create_from_wkt fargs′[3] = kw(:out_warnings, :C_NULL) fargs′[4] = kw(:out_grammar_errors, :C_NULL) diff --git a/src/Proj.jl b/src/Proj.jl index 5b0c69e..7a8fedf 100644 --- a/src/Proj.jl +++ b/src/Proj.jl @@ -4,6 +4,7 @@ using PROJ_jll using CEnum using CoordinateTransformations using NetworkOptions: ca_roots +import GeoFormatTypes as GFT export PROJ_jll export PJ_DIRECTION, PJ_FWD, PJ_IDENT, PJ_INV @@ -88,6 +89,7 @@ const PROJ_COMPUTE_VERSION = VersionNumber const GEODESIC_VERSION_NUM = VersionNumber include("libproj.jl") +include("crs.jl") include("coord.jl") include("error.jl") diff --git a/src/coord.jl b/src/coord.jl index 4c0997a..913f100 100644 --- a/src/coord.jl +++ b/src/coord.jl @@ -115,6 +115,17 @@ function Transformation( return Transformation(pj, direction) end +function Transformation( + source_crs::CRS, + target_crs::CRS; + always_xy::Bool = false, + direction::PJ_DIRECTION = PJ_FWD, + area::Ptr{PJ_AREA} = C_NULL, + ctx::Ptr{PJ_CONTEXT} = C_NULL, +) + return Transformation(source_crs.pj, target_crs.pj; always_xy, direction, area, ctx) +end + function Base.show(io::IO, trans::Transformation) @@ -225,12 +236,31 @@ rectangle. See [`proj_trans_bounds`](https://proj.org/development/reference/functions.html#c.proj_trans_bounds) """ -function bounds(trans::Transformation, (xmin, xmax), (ymin, ymax); densify_pts=21, ctx::Ptr{PJ_CONTEXT}=C_NULL) +function bounds( + trans::Transformation, + (xmin, xmax), + (ymin, ymax); + densify_pts = 21, + ctx::Ptr{PJ_CONTEXT} = C_NULL, +) out_xmin = Ref{Float64}(NaN) out_xmax = Ref{Float64}(NaN) out_ymin = Ref{Float64}(NaN) out_ymax = Ref{Float64}(NaN) - proj_trans_bounds(ctx, trans.pj, trans.direction, xmin, ymin, xmax, ymax, out_xmin, out_ymin, out_xmax, out_ymax, densify_pts) + proj_trans_bounds( + ctx, + trans.pj, + trans.direction, + xmin, + ymin, + xmax, + ymax, + out_xmin, + out_ymin, + out_xmax, + out_ymax, + densify_pts, + ) return (out_xmin[], out_xmax[]), (out_ymin[], out_ymax[]) end diff --git a/src/crs.jl b/src/crs.jl new file mode 100644 index 0000000..c0b7869 --- /dev/null +++ b/src/crs.jl @@ -0,0 +1,139 @@ +""" + CRS(crs) + +Create a CRS. `crs` can be: +- a proj-string, +- a WKT string, +- an object code (like "EPSG:4326", "urn:ogc:def:crs:EPSG::4326", "urn:ogc:def:coordinateOperation:EPSG::1671"), +- an Object name. e.g "WGS 84", "WGS 84 / UTM zone 31N". In that case as uniqueness is not guaranteed, heuristics are applied to determine the appropriate best match. +- a OGC URN combining references for compound coordinate reference systems (e.g "urn:ogc:def:crs,crs:EPSG::2393,crs:EPSG::5717" or custom abbreviated syntax "EPSG:2393+5717"), +- a OGC URN combining references for concatenated operations (e.g. "urn:ogc:def:coordinateOperation,coordinateOperation:EPSG::3895,coordinateOperation:EPSG::1618") +- a PROJJSON string. The jsonschema is at https://proj.org/schemas/v0.4/projjson.schema.json +- a compound CRS made from two object names separated with " + ". e.g. "WGS 84 + EGM96 height" +- a GeoFormatTypes CoordinateReferenceSystemFormat such as EPSG or ProjString +""" +mutable struct CRS + pj::Ptr{PJ} + function CRS(crs::Ptr{PJ}) + crs = new(crs) + finalizer(crs) do crs + crs = proj_destroy(crs) + end + return crs + end +end + +function CRS(crs::AbstractString, ctx::Ptr{PJ_CONTEXT} = C_NULL) + crs = proj_create(crs, ctx) + @assert Bool(proj_is_crs(crs)) "Not a CRS:\n$crs" + return CRS(crs) +end + +function CRS(crs::GFT.CoordinateReferenceSystemFormat, ctx::Ptr{PJ_CONTEXT} = C_NULL) + crs = proj_create(convert(String, crs), ctx) + return CRS(crs) +end + +function Base.show(io::IO, crs::CRS) + info = proj_pj_info(crs) + description = unsafe_string(info.description) + definition = unsafe_string(info.definition) + print( + io, + """CRS + description: $description + definition: $definition + """, + ) +end + +Base.unsafe_convert(::Type{Ptr{Cvoid}}, c::CRS) = c.pj + +function is_type(crs::CRS, types::NTuple{N,PJ_TYPE}) where {N} + if is_compound(crs) + mapreduce(Fix2(is_type, types), |, crs, init = false) + elseif is_bound(crs) + is_type(proj_get_source_crs(crs), types) + else + proj_get_type(crs) in types + end +end + +function Base.iterate(crs::CRS, i = 0) + is_compound(crs) || return nothing + pt = proj_crs_get_sub_crs(crs, i) + if pt == C_NULL + return nothing + else + return CRS(pt), i + 1 + end +end +Base.IteratorSize(::Type{CRS}) = Base.SizeUnknown() +Base.eltype(::Type{CRS}) = CRS + +function is_geographic(crs::CRS) + is_type( + crs, + (PJ_TYPE_GEOGRAPHIC_CRS, PJ_TYPE_GEOGRAPHIC_2D_CRS, PJ_TYPE_GEOGRAPHIC_3D_CRS), + ) +end + +function is_projected(crs::CRS) + is_type(crs, (PJ_TYPE_PROJECTED_CRS,)) +end + +function is_compound(crs::CRS) + proj_get_type(crs) == PJ_TYPE_COMPOUND_CRS +end + +function is_bound(crs::CRS) + proj_get_type(crs) == PJ_TYPE_BOUND_CRS +end + +function GFT.WellKnownText2( + crs::CRS; + type::PJ_WKT_TYPE = PJ_WKT2_2019, + ctx::Ptr{PJ_CONTEXT} = C_NULL, +) + return GFT.WellKnownText2(GFT.CRS(), proj_as_wkt(crs, type, ctx)) +end + +function GFT.WellKnownText( + crs::CRS; + type::PJ_WKT_TYPE = PJ_WKT1_GDAL, + ctx::Ptr{PJ_CONTEXT} = C_NULL, +) + return GFT.WellKnownText(GFT.CRS(), proj_as_wkt(crs, type, ctx)) +end + +function GFT.ESRIWellKnownText( + crs::CRS; + type::PJ_WKT_TYPE = PJ_WKT1_ESRI, + ctx::Ptr{PJ_CONTEXT} = C_NULL, +) + return GFT.ESRIWellKnownText(GFT.CRS(), proj_as_wkt(crs, type, ctx)) +end + +function GFT.ProjString( + crs::CRS; + type::PJ_PROJ_STRING_TYPE = PJ_PROJ_5, + ctx::Ptr{PJ_CONTEXT} = C_NULL, +) + return GFT.ProjString(proj_as_proj_string(crs, type, ctx)) +end + +function GFT.ProjJSON(crs::CRS; ctx::Ptr{PJ_CONTEXT} = C_NULL) + return GFT.ProjJSON(proj_as_projjson(crs, ctx)) +end + +function GFT.EPSG(crs::CRS) + str = proj_get_id_code(crs) + code = parse(Int, str) + return GFT.EPSG(code) +end + +Base.convert(T::Type{<:GFT.CoordinateReferenceSystemFormat}, crs::CRS) = T(crs) +Base.convert(::Type{CRS}, crs::GFT.CoordinateReferenceSystemFormat) = CRS(crs) + +# Maybe enable later, based on https://github.com/JuliaGeo/GeoFormatTypes.jl/issues/21 +# Base.convert(T::Type{<:GFT.CoordinateReferenceSystemFormat}, crs::GFT.CoordinateReferenceSystemFormat) = T(CRS(crs)) diff --git a/src/libproj.jl b/src/libproj.jl index 970ea57..11caf55 100644 --- a/src/libproj.jl +++ b/src/libproj.jl @@ -125,6 +125,11 @@ end """Callback to resolve a filename to a full path""" const proj_file_finder = Ptr{Cvoid} +""" + proj_context_set_file_finder(finder, user_data, ctx = C_NULL) + +` ` +""" function proj_context_set_file_finder(finder, user_data, ctx = C_NULL) @ccall libproj.proj_context_set_file_finder( ctx::Ptr{PJ_CONTEXT}, @@ -151,7 +156,7 @@ end """ proj_context_use_proj4_init_rules(enable, ctx = C_NULL) -Doxygen\\_Suppress +` Doxygen_Suppress ` """ function proj_context_use_proj4_init_rules(enable, ctx = C_NULL) @ccall libproj.proj_context_use_proj4_init_rules( @@ -371,7 +376,7 @@ end """ proj_create(definition, ctx = C_NULL) -Doxygen\\_Suppress +` Doxygen_Suppress ` """ function proj_create(definition, ctx = C_NULL) @ccall libproj.proj_create(ctx::Ptr{PJ_CONTEXT}, definition::Cstring)::Ptr{PJ} @@ -410,6 +415,11 @@ function proj_create_crs_to_crs_from_pj( )::Ptr{PJ} end +""" + proj_normalize_for_visualization(obj, ctx = C_NULL) + +` ` +""" function proj_normalize_for_visualization(obj, ctx = C_NULL) @ccall libproj.proj_normalize_for_visualization( ctx::Ptr{PJ_CONTEXT}, @@ -420,7 +430,7 @@ end """ proj_assign_context(pj, ctx) -Doxygen\\_Suppress +` Doxygen_Suppress ` """ function proj_assign_context(pj, ctx) @ccall libproj.proj_assign_context(pj::Ptr{PJ}, ctx::Ptr{PJ_CONTEXT})::Cvoid @@ -516,6 +526,11 @@ function proj_trans_generic(P, direction, x, sx, nx, y, sy, ny, z, sz, nz, t, st )::Csize_t end +""" + proj_trans_bounds(context, P, direction, xmin, ymin, xmax, ymax, out_xmin, out_ymin, out_xmax, out_ymax, densify_pts) + +` ` +""" function proj_trans_bounds( context, P, @@ -549,7 +564,7 @@ end """ proj_coord(x, y, z = 0.0, t = Inf) -Doxygen\\_Suppress +` Doxygen_Suppress ` """ function proj_coord(x, y, z = 0.0, t = Inf) @ccall libproj.proj_coord(x::Cdouble, y::Cdouble, z::Cdouble, t::Cdouble)::Coord @@ -1056,6 +1071,11 @@ end const PJ_OBJ_LIST = Cvoid +""" + proj_string_list_destroy(list) + +` ` +""" function proj_string_list_destroy(list) @ccall libproj.proj_string_list_destroy(list::Ptr{Cstring})::Cvoid end @@ -1245,11 +1265,11 @@ function proj_get_name(obj) aftercare(@ccall(libproj.proj_get_name(obj::Ptr{PJ})::Cstring)) end -function proj_get_id_auth_name(obj, index) +function proj_get_id_auth_name(obj, index = 0) aftercare(@ccall(libproj.proj_get_id_auth_name(obj::Ptr{PJ}, index::Cint)::Cstring)) end -function proj_get_id_code(obj, index) +function proj_get_id_code(obj, index = 0) aftercare(@ccall(libproj.proj_get_id_code(obj::Ptr{PJ}, index::Cint)::Cstring)) end @@ -1436,6 +1456,11 @@ end const PJ_INSERT_SESSION = Cvoid +""" + proj_insert_object_session_create(ctx = C_NULL) + +` ` +""" function proj_insert_object_session_create(ctx = C_NULL) @ccall libproj.proj_insert_object_session_create( ctx::Ptr{PJ_CONTEXT}, @@ -1497,6 +1522,11 @@ end const PJ_OPERATION_FACTORY_CONTEXT = Cvoid +""" + proj_create_operation_factory_context(authority, ctx = C_NULL) + +` ` +""" function proj_create_operation_factory_context(authority, ctx = C_NULL) @ccall libproj.proj_create_operation_factory_context( ctx::Ptr{PJ_CONTEXT}, @@ -2003,10 +2033,10 @@ end The struct containing information about the ellipsoid. This must be initialized by [`geod_init`](@ref)() before use.******************************************************************** -| Field | Note | -| :---- | :--------------------- | -| a | the equatorial radius | -| f | the flattening SKIP | +| Field | Note | +| :---- | :------------------------ | +| a | the equatorial radius | +| f | the flattening ` SKIP ` | """ struct geod_geodesic a::Cdouble @@ -2028,18 +2058,18 @@ end The struct containing information about a single geodesic. This must be initialized by [`geod_lineinit`](@ref)(), [`geod_directline`](@ref)(), [`geod_gendirectline`](@ref)(), or [`geod_inverseline`](@ref)() before use.******************************************************************** -| Field | Note | -| :---- | :--------------------------------- | -| lat1 | the starting latitude | -| lon1 | the starting longitude | -| azi1 | the starting azimuth | -| a | the equatorial radius | -| f | the flattening | -| salp1 | sine of *azi1* | -| calp1 | cosine of *azi1* | -| a13 | arc length to reference point | -| s13 | distance to reference point SKIP | -| caps | the capabilities | +| Field | Note | +| :---- | :------------------------------------- | +| lat1 | the starting latitude | +| lon1 | the starting longitude | +| azi1 | the starting azimuth | +| a | the equatorial radius | +| f | the flattening | +| salp1 | sine of *azi1* | +| calp1 | cosine of *azi1* | +| a13 | arc length to reference point | +| s13 | distance to reference point ` SKIP ` | +| caps | the capabilities | """ struct geod_geodesicline lat1::Cdouble @@ -2085,11 +2115,11 @@ end The struct for accumulating information about a geodesic polygon. This is used for computing the perimeter and area of a polygon. This must be initialized by [`geod_polygon_init`](@ref)() before use.******************************************************************** -| Field | Note | -| :---- | :--------------------------- | -| lat | the current latitude | -| lon | the current longitude SKIP | -| num | the number of points so far | +| Field | Note | +| :---- | :------------------------------- | +| lat | the current latitude | +| lon | the current longitude ` SKIP ` | +| num | the number of points so far | """ struct geod_polygon lat::Cdouble diff --git a/test/libproj.jl b/test/libproj.jl index c5da2ad..9fd5c40 100644 --- a/test/libproj.jl +++ b/test/libproj.jl @@ -2,6 +2,7 @@ using Test using StaticArrays using Proj import PROJ_jll +import GeoFormatTypes as GFT @testset "Error handling" begin @test Proj.proj_errno_string(0) === nothing @@ -121,7 +122,7 @@ end trans = Proj.Transformation( "EPSG:4326", "+proj=utm +zone=32 +datum=WGS84", - always_xy=true, + always_xy = true, ) # for custom / proj strings, or modified axis order, no description can be looked up in @@ -133,8 +134,8 @@ end direction: forward """ - trans⁻¹ = inv(trans, always_xy=true) - trans¹ = inv(trans⁻¹, always_xy=true) + trans⁻¹ = inv(trans, always_xy = true) + trans¹ = inv(trans⁻¹, always_xy = true) # inv does not flip source and target, so the WKT stays the same wkt_type = Proj.PJ_WKT2_2019 @@ -153,10 +154,10 @@ end trans = Proj.Transformation( "EPSG:4326", "+proj=utm +zone=32 +datum=WGS84", - direction=PJ_IDENT, + direction = PJ_IDENT, ) @test trans(a) === a - trans⁻¹ = inv(trans, always_xy=true) + trans⁻¹ = inv(trans, always_xy = true) @test trans⁻¹(a) === a @test trans⁻¹.direction == PJ_IDENT trans⁻¹.direction = PJ_FWD @@ -187,7 +188,7 @@ end @test is_approx(b, (155191.3538124342, 463537.1362732911)) # with always_xy = true, we need to use lon/lat, and still get x/y out - trans = Proj.Transformation(source_crs, target_crs, always_xy=true) + trans = Proj.Transformation(source_crs, target_crs, always_xy = true) a = (5.39, 52.16) b = Proj.proj_trans(trans.pj, Proj.PJ_FWD, a) @test is_approx(b, (155191.3538124342, 463537.1362732911)) @@ -204,7 +205,8 @@ end @test txmin < tx < txmax @test tymin < ty < tymax - (ttxmin, ttxmax), (ttymin, ttymax) = Proj.bounds(inv(trans), (txmin, txmax), (tymin, tymax)) + (ttxmin, ttxmax), (ttymin, ttymax) = + Proj.bounds(inv(trans), (txmin, txmax), (tymin, tymax)) @test ttxmin < xmin < xmax < ttxmax @test ttymin < ymin < ymax < ttymax @@ -213,7 +215,7 @@ end @testset "dense 4D coord vector transformation" begin source_crs = Proj.proj_create("EPSG:4326") target_crs = Proj.proj_create("EPSG:28992") - trans = Proj.Transformation(source_crs, target_crs, always_xy=true) + trans = Proj.Transformation(source_crs, target_crs, always_xy = true) # This array is mutated in place. Note that this array needs to have 4D elements, # with 2D elements it will only do every other one A = [Proj.Coord(5.39, 52.16) for _ = 1:5] @@ -235,7 +237,7 @@ end @testset "generic array transformation" begin source_crs = Proj.proj_create("EPSG:4326") target_crs = Proj.proj_create("EPSG:28992") - trans = Proj.Transformation(source_crs, target_crs, always_xy=true) + trans = Proj.Transformation(source_crs, target_crs, always_xy = true) # inplace transformation of vector of 2D coordinates # using https://proj.org/development/reference/functions.html#c.proj_trans_generic @@ -266,8 +268,8 @@ end end @testset "compose" begin - trans1 = Proj.Transformation("EPSG:4326", "EPSG:28992", always_xy=true) - trans2 = Proj.Transformation("EPSG:32632", "EPSG:2027", always_xy=true) + trans1 = Proj.Transformation("EPSG:4326", "EPSG:28992", always_xy = true) + trans2 = Proj.Transformation("EPSG:32632", "EPSG:2027", always_xy = true) trans = trans1 ∘ trans2 # same as compose(trans1, trans2) source_crs = Proj.proj_get_source_crs(trans.pj) target_crs = Proj.proj_get_target_crs(trans.pj) @@ -320,7 +322,7 @@ end end @testset "in and output types" begin - trans = Proj.Transformation("EPSG:4326", "EPSG:28992", always_xy=true) + trans = Proj.Transformation("EPSG:4326", "EPSG:28992", always_xy = true) trans(Proj.proj_coord(5.39, 52.16)) b = trans(SA[5.39, 52.16, 0.0, 0.0]) @@ -351,12 +353,12 @@ end # turn off network, no z transformation @test !Proj.enable_network!(false) @test !Proj.network_enabled() - trans_z = Proj.Transformation("EPSG:4326+5773", "EPSG:7856+5711", always_xy=true) + trans_z = Proj.Transformation("EPSG:4326+5773", "EPSG:7856+5711", always_xy = true) @test trans_z((151, -33, 5))[3] == 5 # turn on network, z transformation @test Proj.enable_network!(true) @test Proj.network_enabled() - trans_z = Proj.Transformation("EPSG:4326+5773", "EPSG:7856+5711", always_xy=true) + trans_z = Proj.Transformation("EPSG:4326+5773", "EPSG:7856+5711", always_xy = true) z = trans_z((151, -33, 5))[3] @test z ≈ 5.280647277836724f0 @@ -368,3 +370,71 @@ end # restore setting as outside the testset Proj.enable_network!(as_before) end + +@testset "CRS" begin + gftcrs = GFT.EPSG(4326) + + crs = Proj.CRS("EPSG:4326") + crs = Proj.CRS("+proj=longlat +datum=WGS84 +no_defs +type=crs") + @test_throws AssertionError Proj.CRS( + "+proj=pipeline +ellps=GRS80 +step +proj=merc +step +proj=axisswap +order=2,1", + ) + + crs = Proj.CRS(gftcrs) + + @test repr(crs) == """ + CRS + description: WGS 84 + definition: \n""" + + @test Proj.is_geographic(crs) + @test !Proj.is_projected(crs) + @test !Proj.is_compound(crs) + @test !Proj.is_bound(crs) + + @test isnothing(iterate(crs)) + + + + trans = Proj.Transformation(crs, crs) + @test repr(trans) == """ + Transformation noop + description: Null geographic offset from WGS 84 to WGS 84 + definition: proj=noop ellps=GRS80 + direction: forward + """ + + crs = Proj.CRS("EPSG:4326+3855") + @test length(collect(crs)) == 2 +end + +@testset "GFT" begin + gftcrs = GFT.EPSG(4326) + crs = Proj.CRS(gftcrs) + + @test GFT.WellKnownText(crs) == GFT.WellKnownText{GFT.CRS}( + GFT.CRS(), + "GEOGCS[\"WGS 84\",\n DATUM[\"WGS_1984\",\n SPHEROID[\"WGS 84\",6378137,298.257223563,\n AUTHORITY[\"EPSG\",\"7030\"]],\n AUTHORITY[\"EPSG\",\"6326\"]],\n PRIMEM[\"Greenwich\",0,\n AUTHORITY[\"EPSG\",\"8901\"]],\n UNIT[\"degree\",0.0174532925199433,\n AUTHORITY[\"EPSG\",\"9122\"]],\n AUTHORITY[\"EPSG\",\"4326\"]]", + ) + @test GFT.WellKnownText2(crs) == GFT.WellKnownText2{GFT.CRS}( + GFT.CRS(), + "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]", + ) + @test GFT.ESRIWellKnownText(crs) == GFT.ESRIWellKnownText{GFT.CRS}( + GFT.CRS(), + "GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",6378137.0,298.257223563]],PRIMEM[\"Greenwich\",0.0],UNIT[\"Degree\",0.0174532925199433]]", + ) + @test GFT.ProjString(crs) == + GFT.ProjString("+proj=longlat +datum=WGS84 +no_defs +type=crs") + @test GFT.ProjJSON(crs) == GFT.ProjJSON( + "{\n \"\$schema\": \"https://proj.org/schemas/v0.5/projjson.schema.json\",\n \"type\": \"GeographicCRS\",\n \"name\": \"WGS 84\",\n \"datum_ensemble\": {\n \"name\": \"World Geodetic System 1984 ensemble\",\n \"members\": [\n {\n \"name\": \"World Geodetic System 1984 (Transit)\",\n \"id\": {\n \"authority\": \"EPSG\",\n \"code\": 1166\n }\n },\n {\n \"name\": \"World Geodetic System 1984 (G730)\",\n \"id\": {\n \"authority\": \"EPSG\",\n \"code\": 1152\n }\n },\n {\n \"name\": \"World Geodetic System 1984 (G873)\",\n \"id\": {\n \"authority\": \"EPSG\",\n \"code\": 1153\n }\n },\n {\n \"name\": \"World Geodetic System 1984 (G1150)\",\n \"id\": {\n \"authority\": \"EPSG\",\n \"code\": 1154\n }\n },\n {\n \"name\": \"World Geodetic System 1984 (G1674)\",\n \"id\": {\n \"authority\": \"EPSG\",\n \"code\": 1155\n }\n },\n {\n \"name\": \"World Geodetic System 1984 (G1762)\",\n \"id\": {\n \"authority\": \"EPSG\",\n \"code\": 1156\n }\n },\n {\n \"name\": \"World Geodetic System 1984 (G2139)\",\n \"id\": {\n \"authority\": \"EPSG\",\n \"code\": 1309\n }\n }\n ],\n \"ellipsoid\": {\n \"name\": \"WGS 84\",\n \"semi_major_axis\": 6378137,\n \"inverse_flattening\": 298.257223563\n },\n \"accuracy\": \"2.0\",\n \"id\": {\n \"authority\": \"EPSG\",\n \"code\": 6326\n }\n },\n \"coordinate_system\": {\n \"subtype\": \"ellipsoidal\",\n \"axis\": [\n {\n \"name\": \"Geodetic latitude\",\n \"abbreviation\": \"Lat\",\n \"direction\": \"north\",\n \"unit\": \"degree\"\n },\n {\n \"name\": \"Geodetic longitude\",\n \"abbreviation\": \"Lon\",\n \"direction\": \"east\",\n \"unit\": \"degree\"\n }\n ]\n },\n \"scope\": \"Horizontal component of 3D system.\",\n \"area\": \"World.\",\n \"bbox\": {\n \"south_latitude\": -90,\n \"west_longitude\": -180,\n \"north_latitude\": 90,\n \"east_longitude\": 180\n },\n \"id\": {\n \"authority\": \"EPSG\",\n \"code\": 4326\n }\n}", + ) + @test GFT.EPSG(crs) == GFT.EPSG("EPSG:4326") + + @test convert(GFT.EPSG, crs) == GFT.EPSG("EPSG:4326") + @test Proj.proj_get_id_code(convert(Proj.CRS, GFT.EPSG("EPSG:4326"))) == + Proj.proj_get_id_code(crs) + + # Maybe enable later, based on https://github.com/JuliaGeo/GeoFormatTypes.jl/issues/21 + # @test convert(GFT.ProjString, gftcrs) == GFT.ProjString("+proj=longlat +datum=WGS84 +no_defs +type=crs") +end