From fccecd35a33b3532ac663da3485ff381ea14a99c Mon Sep 17 00:00:00 2001 From: rafaqz Date: Fri, 5 May 2023 16:03:10 +0200 Subject: [PATCH] move to transformations --- src/GeometryOps.jl | 3 +- src/transformations/flip.jl | 18 ++++++++++ src/transformations/reproject.jl | 57 ++++++++++++++++++++++++++++++++ test/methods/reproject.jl | 39 ---------------------- test/transformations.jl | 52 +++++++++++++++++++++++++++++ 5 files changed, 129 insertions(+), 40 deletions(-) create mode 100644 src/transformations/flip.jl create mode 100644 src/transformations/reproject.jl delete mode 100644 test/methods/reproject.jl create mode 100644 test/transformations.jl diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index c0f20dc91..93839d47e 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -11,6 +11,7 @@ include("methods/signed_distance.jl") include("methods/signed_area.jl") include("methods/centroid.jl") include("methods/contains.jl") -include("methods/reproject.jl") +include("transformations/reproject.jl") +include("transformations/flip.jl") end diff --git a/src/transformations/flip.jl b/src/transformations/flip.jl new file mode 100644 index 000000000..cdbedc8fa --- /dev/null +++ b/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 diff --git a/src/transformations/reproject.jl b/src/transformations/reproject.jl new file mode 100644 index 000000000..7986c21f2 --- /dev/null +++ b/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 diff --git a/test/methods/reproject.jl b/test/methods/reproject.jl deleted file mode 100644 index ee2324c9f..000000000 --- a/test/methods/reproject.jl +++ /dev/null @@ -1,39 +0,0 @@ -using Test - -import GeoInterface as GI -import GeometryOps as GO -using GeoFormatTypes -import Proj - -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) diff --git a/test/transformations.jl b/test/transformations.jl new file mode 100644 index 000000000..3062c302c --- /dev/null +++ b/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