add geointerface reprojection

Project.toml

@@ -6,13 +6,15 @@ version = "1.4.0"
 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"
 PROJ_jll = "900.100"
 julia = "1.6"
 

src/Proj.jl

@@ -3,6 +3,7 @@ module Proj
 using PROJ_jll
 using CEnum
 using CoordinateTransformations
+using GeoInterface
 using NetworkOptions: ca_roots
 import GeoFormatTypes as GFT
 
@@ -92,6 +93,7 @@ include("libproj.jl")
 include("crs.jl")
 include("coord.jl")
 include("error.jl")
+include("geointerface.jl")
 
 """
     unsafe_loadstringlist(ptr::Ptr{Cstring})

src/geointerface.jl

@@ -0,0 +1,98 @@
+using GeoInterface
+
+const GI = GeoInterface
+
+# Coord is a PointTrait
+GI.isgeometry(::Type{Coord}) = true
+GI.geomtrait(::Coord) = GI.PointTrait()
+GI.x(::PointTrait, c::Coord) = c.x
+GI.y(::PointTrait, c::Coord) = c.y
+GI.z(::PointTrait, c::Coord) = c.z
+GI.getcoord(::PointTrait, c::Coord) = (c.x, c.y, c.z)
+GI.getcoord(::PointTrait, c::Coord, i::Int) = c[i]
+coordnames(::PointTrait, ::Coord) = (:X, :Y, :Z)
+
+"""
+    reproject(geometry, source_crs, target_crs)
+
+Reproject any GeoInterface.jl compatible `geometry` from `source_crs` to `target_crs`.
+
+The returned object will be constructed from `GeoInterface.WrapperGeometry`
+geometries, wrapping views of `Proj.Coord`.
+"""
+function reproject(geom, source_crs, target_crs; time=Inf)
+    source_crs1 = convert(Proj.CRS, source_crs)
+    target_crs1 = convert(Proj.CRS, target_crs)
+    trans = Proj.Transformation(source_crs1, target_crs1, always_xy = true)
+    coords = if GI.is3d(geom)
+        [Proj.Coord(GI.x(p), GI.y(p), GI.z(p), time) for p in GI.getpoint(geom)]
+    else
+        [Proj.Coord(GI.x(p), GI.y(p), 0.0, time) for p in GI.getpoint(geom)]
+    end
+    err = Proj.proj_trans_array(trans.pj, Proj.PJ_FWD, length(coords), coords)
+    err == 0 || error("Proj error $err")
+
+    crs = target_crs isa GFT.GeoFormat ? target_crs : convert(WellKnownText, crs)
+    return reconstruct(geom, coords; crs)
+end
+
+
+"""
+    reconstruct(geometry, points::AbstractVector)
+
+Reconstruct GeoInterface compatible `geometry` from the `PointTrait` geometries in `points`.
+
+The returned object will be constructed from `GeoInterface.WrapperGeometry`
+geometries, wrapping views into the `points` vector. Otherwise, the structure and
+GeoInterface traits will be the same as for `geometry`.
+
+This is a lazy operation, only allocating for outer wrappers of nested geometries.
+Later changes to the `points` vector *will* affect the returned geometry.
+
+`GeoInterface.npoint(geometry) == length(points)` must be `true`.
+"""
+function reconstruct(geom, points::AbstractVector; crs=nothing)
+    T = nonmissingtype(eltype(points))
+    isgeometry(T) || throw(ArgumentError("points of type $(T) are not GeoInterface.jl compatible objects"))
+    trait = GI.trait(first(skipmissing(points))) 
+    trait isa PointTrait || throw(ArgumentError("can only reconstruct from an array of points, got $trait"))
+    n, geom = _reconstruct(geom, points, 0, crs)
+    @assert n == GI.npoint(geom)
+    return geom
+end
+
+_reconstruct(geom, points, n::Int, crs) = _reconstruct(GI.trait(geom), geom, points, n::Int, crs)
+# Nested geometries. We need to reconstruct their child geoms.
+function _reconstruct(trait::AbstractGeometryTrait, geom, points, n, crs)
+    T = GI.geointerface_geomtype(trait)
+    geoms = map(GI.getgeom(geom)) do childgeom
+        n, reconstructed_childgeom = _reconstruct(GI.trait(childgeom), childgeom, points, n, crs)
+        reconstructed_childgeom
+    end
+    return n, T(geoms; crs)
+end
+# Bottom-level geometries that can wrap a vector of points.
+function _reconstruct(
+    trait::Union{GI.LineTrait,GI.LineStringTrait,GI.LinearRingTrait,GI.MultiPointTrait}, 
+    geom, points, n, crs
+)
+    T = GI.geointerface_geomtype(trait)
+    npoints = GI.npoint(geom)
+    v = view(points, n + 1:n + npoints)
+    geom = if GI.is3d(geom)
+        GI.ismeasured(geom) ? T{true,true}(v; crs) : T{true,false}(v; crs)
+    else
+        GI.ismeasured(geom) ? T{false,true}(v; crs) : T{false,false}(v; crs)
+    end
+    return n + npoints, geom
+end
+# Points
+function _reconstruct(::GI.PointTrait, geom, points, n, crs)
+    n += 1
+    point = if GI.is3d(geom)
+        GI.ismeasured(geom) ? GI.Point{true,true}(points[n]) : GI.Point{true,false}(points[n])
+    else
+        GI.ismeasured(geom) ? GI.Point{false,true}(points[n]) : GI.Point{false,false}(points[n])
+    end
+    return n, point
+end

test/geointerface.jl

@@ -0,0 +1,14 @@
+using Proj
+using GeoInterface
+using GeoFormatTypes
+const GI = GeoInterface
+
+linearring1 = GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)])
+linearring2 = GI.LinearRing([(11, 2), (13, 4), (15, 6), (11, 2)])
+
+polygon1 = GI.Polygon([linearring1])
+polygon2 = GI.Polygon([linearring2])
+multipolygon = GI.MultiPolygon([polygon1, polygon2])
+
+Proj.reproject(multipolygon, EPSG(4326), EPSG(3857))
+