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,172 @@
+using GeoInterface
+
+const GI = GeoInterface
+
+struct Point{D}
+    coord::Coord
+end
+Point{D}(x, y, z, t) where D = Point{D}(Coord(x, y, z, y))
+
+Point{3}(x, y, z=0.0; time=Inf) = Point{3}(x, y, z, time)
+Point{2}(x, y; time=Inf) = Point{2}(x, y, 0.0, time)
+Point(x, y; time=Inf) = Point{2}(x, y, 0.0, time)
+Point(x, y, z; time=Inf) = Point{3}(x, y, z, time)
+Point(x, y, z, t) = Point{3}(x, y, z, t)
+
+# this shields a StackOverflow from the splatting constructor
+Point(::Real) = error("Proj.Point takes 2 to 4 numbers, one given")
+Point(v) = Point(v...)
+function Point(v::AbstractVector{<:Real})
+    n = length(v)
+    return if n == 2
+        Point{2}(v[begin], v[begin+1])
+    elseif n == 3
+        Point{3}(v[begin], v[begin+1], v[begin+2])
+    elseif n == 4
+        Point{3}(v[begin], v[begin+1], v[begin+2], v[begin+3])
+    else
+        error("Proj.Point takes 2 to 4 numbers")
+    end
+end
+
+GI.isgeometry(::Type{<:Point}) = true
+GI.geomtrait(::Point) = GI.PointTrait()
+GI.x(::PointTrait, c::Point) = c.coord.x # TODO: these may have the wrong order for some CRS. 
+GI.y(::PointTrait, c::Point) = c.coord.y 
+GI.z(::PointTrait, c::Point{3}) = c.coord.z
+GI.is3d(::PointTrait, c::Point{2}) = false
+GI.is3d(::PointTrait, c::Point{3}) = true
+GI.ismeasured(::PointTrait, c::Point) = false
+function GI.getcoord(::PointTrait, p::Point{N}, i::Int) where N
+    checkbounds(1:N, i)
+    @inbounds p.coord[i]
+end
+GI.getcoord(::PointTrait, p::Point{2}) = (c.x, c.y)
+GI.getcoord(::PointTrait, p::Point{3}) = (c.x, c.y, c.z)
+coordnames(::PointTrait, ::Point{2}) = (:X, :Y)
+coordnames(::PointTrait, ::Point{3}) = (:X, :Y, :Z)
+
+"""
+    reproject(geometry; source_crs, target_crs, transform, always_xy, time)
+    reproject(geometry, source_crs, target_crs; always_xy, time)
+    reproject(geometry, transform; always_xy, time)
+
+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 a `Vector{Proj.Point{D}}`, where `D` is the dimension.
+
+# Arguments
+
+- `geometry`: Any GeoInterface.jl compatible geometries.
+- `source_crs`: the source coordinate referece system, as a GeoFormatTypes.jl object or a string.
+- `target_crs`: the target coordinate referece system, as a GeoFormatTypes.jl object or a string.
+
+If these a passed as keywords, `transform` will take priority. 
+Without it `target_crs` is always needed, and `source_crs` is
+needed if it is not retreivable from the geometry with `GeoInterface.crs(geometry)`.
+
+# Keywords
+
+-`always_xy`: force x, y coordinate order, `true` by default. 
+    `false` will expect and return points in the crs coordinate order.
+-`time`: the time for the coordinates. `Inf` by default.
+"""
+function reproject(geom; 
+    source_crs=nothing, target_crs=nothing, transform=nothing, kw...
+)
+    if isnothing(transform)
+        source_crs = isnothing(source_crs) ? GeoInterface.crs(geom) : source_crs
+        isnothing(source_crs) && throw(ArgumentError("geom has no crs attatched. Pass a `source_crs` keyword"))
+        reproject(geom, source_crs, target_crs; transform, kw...)
+    else
+        reproject(geom, transform)
+    end
+end
+function reproject(geom, source_crs, target_crs; 
+    time=Inf, 
+    always_xy=true,
+    transform=Proj.Transformation(source_crs1, target_crs1; always_xy),
+)
+    reproject(geom, transform; time)
+end
+function reproject(geom, transform::Transformation; time)
+    wrapped_geom, coords = if GI.is3d(geom)
+        reconstruct(geom; crs=target_crs) do p
+            Proj.Point{3}(GI.x(p), GI.y(p), GI.z(p); time)
+        end
+    else
+        reconstruct(geom; crs=target_crs) do p
+            Proj.Point{2}(GI.x(p), GI.y(p); time)
+        end
+    end
+    source_crs1 = convert(Proj.CRS, source_crs)
+    target_crs1 = convert(Proj.CRS, target_crs)
+    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 wrapped_geom
+end
+
+
+    # reconstruct(f::Function, geom)
+
+# Reconstruct GeoInterface compatible `geometry` from the `PointTrait` geometries
+# in `points`, or from the result of function `f` over the points in `geom`.
+
+# 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.
+function reconstruct(f, geom; crs=nothing)
+    p1 = f(first(GI.getpoint(geom)))
+    T = typeof(p1)
+    isgeometry(T) || throw(ArgumentError("points of type $(T) are not GeoInterface.jl compatible objects"))
+    trait = GI.trait(p1)
+    trait isa PointTrait || throw(ArgumentError("can only reconstruct from an array of points, got $trait"))
+    points = Vector{T}(undef, GI.npoint(geom))
+    n, geom = _reconstruct!(f, points, geom, 0, crs)
+    @assert n == GI.npoint(geom)
+    return geom, points
+end
+
+_reconstruct!(f, points, geom, n::Int, crs) = _reconstruct!(f, points, GI.trait(geom), geom, n::Int, crs)
+# Nested geometries. We need to reconstruct their child geoms.
+function _reconstruct!(f, points, trait::AbstractGeometryTrait, geom, n, crs)
+    T = GI.geointerface_geomtype(trait)
+    geoms = map(GI.getgeom(geom)) do childgeom
+        childcrs = nothing # We dont pass the CRS down, the type gets too complicated
+        n, reconstructed_childgeom = _reconstruct!(f, points, GI.trait(childgeom), childgeom, n, crs)
+        reconstructed_childgeom
+    end
+    return n, T(geoms; crs)
+end
+# Bottom-level geometries that can wrap a vector of points.
+function _reconstruct!(f, points, 
+    trait::Union{GI.LineTrait,GI.LineStringTrait,GI.LinearRingTrait,GI.MultiPointTrait}, 
+    geom, n, crs
+)
+    T = GI.geointerface_geomtype(trait)
+    npoints = GI.npoint(geom)
+    v = view(points, n + 1:n + npoints)
+    v .= f.(GI.getpoint(geom))
+    wrapped_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, wrapped_geom
+end
+# Points
+function _reconstruct!(f, points, ::GI.PointTrait, geom, n, crs)
+    n += 1
+    points[n] = f(point)
+    wrapped_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, wrapped_point
+end

test/geointerface.jl

@@ -0,0 +1,40 @@
+using Proj
+using Test
+using GeoInterface
+using GeoFormatTypes
+const GI = GeoInterface
+
+using GeoInterfaceRecipes
+GeoInterfaceRecipes.@enable_geo_plots GeoInterface.Wrappers.WrapperGeometry
+
+ring1 = GI.LinearRing([(1, 2), (7, 4), (5, 6), (1, 2)])
+ring2 = GI.LinearRing([(11, 2), (20, 4), (15, 6), (11, 2)])
+hole2 = GI.LinearRing([(14, 4), (16, 4), (17, 5), (14, 4)])
+
+# Set up a regular tranformation of the points for reference
+source_crs = convert(Proj.CRS, EPSG(4326))
+target_crs = convert(Proj.CRS, EPSG(3857))
+trans = Proj.Transformation(source_crs, target_crs; always_xy=true)
+ref_points3857 = map(GeoInterface.getpoint(multipolygon)) do p
+    trans([GeoInterface.x(p), GeoInterface.y(p)])
+end
+
+polygon1 = GI.Polygon([ring1])
+polygon2 = GI.Polygon([ring2, hole2])
+multipolygon = GI.MultiPolygon([polygon1, polygon2])
+
+multipolygon3857 = Proj.reproject(multipolygon, EPSG(4326), EPSG(3857))
+multipolygon4326 = Proj.reproject(multipolygon3857; target_crs=EPSG(4326))
+points4326_1 = collect(GI.getpoint(multipolygon))
+points4326_2 = GI.convert.(Tuple, GI.getpoint(multipolygon4326))
+points3857 = GI.convert.(Tuple, GI.getpoint(multipolygon3857))
+
+# Comparison to regular `trans` on points
+@test all(map((p1, p2) -> all(map(isapprox, p1, p2)), ref_points3857, points3857))
+
+# Round trip comparison
+@test all(map((p1, p2) -> all(map(isapprox, p1, p2)), points4326_1, points4326_2))
+
+# Embedded crs check
+@test GI.crs(multipolygon3857) == EPSG(3857)
+@test GI.crs(multipolygon4326) == EPSG(4326)

test/runtests.jl

@@ -8,5 +8,6 @@ is_approx(a, b) = all(isapprox(c[1], c[2]) for c in zip(a, b))
 @testset "Proj" begin
     include("libproj.jl")
     include("applications.jl")
+    include("geointerface.jl")
 
 end # testset "Proj"