Implement GeoInterface for Coord and in Transformations. (#91)

* Implement GeoInterface for Coord and in Transformations.

* Proj.identify (#92)

* add high level interface to proj_identify

* add tests for Proj.identify

* move identify from libproj to crs

* address comments on Proj.identify

* update tests

* free pj_list

* final improvements

- standardize docstring
- avoid intermediate Vectors
- free out_confidence

* run JuliaFormatter

---------

Co-authored-by: Martijn Visser <[email protected]>

* run JuliaFormatter

---------

Co-authored-by: Alex Gardner <[email protected]>
Co-authored-by: Martijn Visser <[email protected]>

Project.toml

@@ -1,18 +1,20 @@
 name = "Proj"
 uuid = "c94c279d-25a6-4763-9509-64d165bea63e"
-version = "1.4.0"
+version = "1.5.0"
 
 [deps]
 CEnum = "fa961155-64e5-5f13-b03f-caf6b980ea82"
 CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
 GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
+GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 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"
+GeoInterface = "1.3"
 PROJ_jll = "900.100"
 julia = "1.6"
 

src/Proj.jl

@@ -5,6 +5,7 @@ using CEnum
 using CoordinateTransformations
 using NetworkOptions: ca_roots
 import GeoFormatTypes as GFT
+import GeoInterface as GI
 
 export PROJ_jll
 export PJ_DIRECTION, PJ_FWD, PJ_IDENT, PJ_INV
@@ -83,6 +84,31 @@ Base.getindex(coord::Coord, i::Int) = getfield(coord, i)
 Base.IndexStyle(::Type{Coord}) = IndexLinear()
 Base.eltype(::Coord) = Float64
 
+geointerface_geomtype(::GI.PointTrait) = Coord
+
+GI.isgeometry(::Type{Coord}) = true
+GI.geomtrait(::Coord) = GI.PointTrait()
+GI.ncoord(::Coord) = 4
+GI.getcoord(coord::Coord, i::Int) = getfield(coord, i)
+GI.x(coord::Coord) = coord.x
+GI.y(coord::Coord) = coord.y
+GI.z(coord::Coord) = coord.z
+GI.m(coord::Coord) = coord.t
+
+function GI.convert(::Type{Coord}, ::GI.PointTrait, geom)
+    n = GI.ncoord(geom)
+    if n == 2
+        return Coord(GI.x(geom), GI.y(geom))
+    elseif n == 3
+        return Coord(GI.x(geom), GI.y(geom), GI.z(geom))
+    elseif n == 4
+        return Coord(GI.x(geom), GI.y(geom), GI.z(geom), GI.m(geom))
+    else
+        error("Coord takes 2 to 4 numbers")
+    end
+end
+
+
 # type aliases
 const NTuple234 = Union{NTuple{2,Float64},NTuple{3,Float64},NTuple{4,Float64}}
 const PROJ_COMPUTE_VERSION = VersionNumber

src/coord.jl

@@ -232,7 +232,7 @@ function (trans::Transformation)(coord::Coord)::NTuple{4,Float64}
     return p.x, p.y, p.z, p.t
 end
 
-function (trans::Transformation)(coord)::NTuple234
+function (trans::Transformation)(coord::NTuple234)::NTuple234
     n = length(coord)
     p = proj_trans(trans.pj, trans.direction, coord)
     return if n == 2
@@ -244,6 +244,21 @@ function (trans::Transformation)(coord)::NTuple234
     end
 end
 
+function (trans::Transformation)(coord)
+    (GI.isgeometry(coord) && GI.geomtrait(coord)) == GI.PointTrait() ||
+        throw(ArgumentError("Argument is not a Point geometry"))
+    c = GI.convert(Coord, GI.PointTrait(), coord)
+    p = proj_trans(trans.pj, trans.direction, c)
+    n = GI.ncoord(coord)
+    if n == 2
+        return p.x, p.y
+    elseif n == 3
+        return p.x, p.y, p.z
+    else
+        return p.x, p.y, p.z, p.t
+    end
+end
+
 """
     Proj.bounds(trans::Transformation, (xmin, xmax), (ymin,ymax); densify_pts=21) -> ((bxmin, bxmax), (bymin, bymax))
 

test/libproj.jl

@@ -1,6 +1,7 @@
 using Test
 using StaticArrays
 using Proj
+import GeoInterface as GI
 import PROJ_jll
 import GeoFormatTypes as GFT
 
@@ -478,6 +479,23 @@ end
     # @test convert(GFT.ProjString, gftcrs) == GFT.ProjString("+proj=longlat +datum=WGS84 +no_defs +type=crs")
 end
 
+@testset "GeoInterface" begin
+    @test GI.testgeometry(Proj.Coord(0, 0))
+
+    wp = GI.Wrappers.Point(5, 54)
+    np = GI.convert(Proj.Coord, wp)
+    @test np isa Proj.Coord
+
+    np = GI.convert(Proj, wp)
+    @test np isa Proj.Coord
+
+    trans = Proj.Transformation("EPSG:4326", "EPSG:28992", always_xy = true)
+    tp = trans(wp)
+    @test length(tp) == 2
+    @test tp[1] ≈ 129604.1711
+    @test tp[2] ≈ 668374.4875
+end
+
 @testset "Proj.identify" begin
     crs = Proj.CRS("+proj=longlat +datum=WGS84 +no_defs +type=crs")
     identity = Proj.identify(crs)