Merge pull request #2 from asinghvi17/reproject

Add a `reproject` and `flip` method

Project.toml

@@ -6,9 +6,11 @@ version = "0.0.1-DEV"
 [deps]
 GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
+Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 
 [compat]
 GeometryBasics = "0.4.7"
+Proj = "1"
 julia = "1.3"
 
 [extras]

src/GeometryOps.jl

@@ -2,6 +2,7 @@ module GeometryOps
 
 using GeoInterface
 using GeometryBasics
+import Proj
 
 const GI = GeoInterface
 
@@ -10,5 +11,7 @@ include("methods/signed_distance.jl")
 include("methods/signed_area.jl")
 include("methods/centroid.jl")
 include("methods/contains.jl")
+include("transformations/reproject.jl")
+include("transformations/flip.jl")
 
 end

src/methods/reproject.jl

@@ -0,0 +1,57 @@
+
+
+"""
+    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; kw...)
+    else
+        reproject(geom, transform; kw...)
+    end
+end
+function reproject(geom, source_crs, target_crs;
+    time=Inf,
+    always_xy=true,
+    transform=Proj.Transformation(Proj.CRS(source_crs), Proj.CRS(target_crs); always_xy),
+)
+    reproject(geom, transform; time, target_crs)
+end
+function reproject(geom, transform::Proj.Transformation; time=Inf, target_crs=nothing)
+    if GI.is3d(geom)
+        return apply(PointTrait, geom; crs=target_crs) do p
+            transform(GI.x(p), GI.y(p), GI.z(p))
+        end
+    else
+        return apply(PointTrait, geom; crs=target_crs) do p
+            transform(GI.x(p), GI.y(p))
+        end
+    end
+end

src/transformations/flip.jl

@@ -0,0 +1,18 @@
+"""
+    flip(obj)
+
+Swap all of the x and y coordinates in obj, otherwise
+keeping the original structure (but not necessarily the
+original type).
+"""
+function flip(geom) 
+    if GI.is3d(geom)
+        return apply(PointTrait, geom) do point
+            (GI.y(p), GI.x(p), GI.z(p))
+        end
+    else
+        return apply(PointTrait, geom) do point
+            (GI.y(p), GI.x(p))
+        end
+    end
+end

src/transformations/reproject.jl

@@ -0,0 +1,57 @@
+
+
+"""
+    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; kw...)
+    else
+        reproject(geom, transform; kw...)
+    end
+end
+function reproject(geom, source_crs, target_crs;
+    time=Inf,
+    always_xy=true,
+    transform=Proj.Transformation(Proj.CRS(source_crs), Proj.CRS(target_crs); always_xy),
+)
+    reproject(geom, transform; time, target_crs)
+end
+function reproject(geom, transform::Proj.Transformation; time=Inf, target_crs=nothing)
+    if GI.is3d(geom)
+        return apply(PointTrait, geom; crs=target_crs) do p
+            transform(GI.x(p), GI.y(p), GI.z(p))
+        end
+    else
+        return apply(PointTrait, geom; crs=target_crs) do p
+            transform(GI.x(p), GI.y(p))
+        end
+    end
+end

test/transformations.jl

@@ -0,0 +1,52 @@
+using Test
+
+import GeoInterface as GI
+import GeometryOps as GO
+using GeoFormatTypes
+import Proj
+
+@testset "reproject" begin
+    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)
+
+    polygon1 = GI.Polygon([ring1])
+    polygon2 = GI.Polygon([ring2, hole2])
+    multipolygon = GI.MultiPolygon([polygon1, polygon2])
+
+    ref_points3857 = map(GI.getpoint(multipolygon)) do p
+        trans([GI.x(p), GI.y(p)])
+    end
+
+    multipolygon3857 = GO.reproject(multipolygon, EPSG(4326), EPSG(3857))
+    multipolygon4326 = GO.reproject(multipolygon3857; target_crs=EPSG(4326))
+    points4326_1 = collect(GI.getpoint(multipolygon))
+    points4326_2 = GI.getcoord.(GI.getpoint(multipolygon4326))
+    points3857 = GI.getcoord.(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)
+end
+
+
+
+@testset "flip" begin
+    geom = GI.Polygon([GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)]), 
+                       GI.LinearRing([(3, 4), (5, 6), (6, 7), (3, 4)])])
+
+
+    @test GO.flip(geom) == GI.Polygon([GI.LinearRing([(2, 1), (4, 3), (6, 5), (2, 1)]), 
+                                       GI.LinearRing([(4, 3), (6, 5), (7, 6), (4, 3)])])
+end