Add CRS constructor.

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"
 

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.

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)

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")
 

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
 

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))