diff --git a/docs/Project.toml b/docs/Project.toml index b18a68afcc..7f3bd6df08 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -42,6 +42,7 @@ NaturalEarth = "436b0209-26ab-4e65-94a9-6526d86fea76" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" Proj = "c94c279d-25a6-4763-9509-64d165bea63e" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc" Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4" SortTileRecursiveTree = "746ee33f-1797-42c2-866d-db2fce69d14d" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" diff --git a/docs/src/api.md b/docs/src/api.md index e5d159cae8..c4c44998c4 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -65,4 +65,5 @@ Modules = [GeometryOps] ```@autodocs Modules = [GeometryOpsCore] +Order = [:type] ``` diff --git a/docs/src/tutorials/creating_geometry.md b/docs/src/tutorials/creating_geometry.md index f248587239..194ef5ef14 100644 --- a/docs/src/tutorials/creating_geometry.md +++ b/docs/src/tutorials/creating_geometry.md @@ -20,6 +20,7 @@ using GeoJSON # to load some data # Packages for coordinate transformation and projection import CoordinateTransformations import Proj +import Rotations # Plotting using CairoMakie using GeoMakie @@ -79,7 +80,7 @@ This time we get a bit more fancy with point creation. ````@example creating_geometry r = 2; k = 10; -ϴ = 0:0.01:2pi; +ϴ = LinRange(0, 2pi, 629); x = r .* (k + 1) .* cos.(ϴ) .- r .* cos.((k + 1) .* ϴ); y = r .* (k + 1) .* sin.(ϴ) .- r .* sin.((k + 1) .* ϴ); lines = GI.LineString(GI.Point.(zip(x,y))); @@ -93,7 +94,8 @@ A `LinearRing` is simply a `LineString` with the same beginning and endpoint, i. A `LinearRing` is composed of a series of points. ````@example creating_geometry -ring1 = GI.LinearRing(GI.getpoint(lines)); +ring1_points = collect(GI.getpoint(lines)) +ring1 = GI.LinearRing(vcat(ring1_points, [first(ring1_points)])); ```` Now, let's make the `LinearRing` into a `Polygon`. @@ -113,6 +115,37 @@ plot!(polygon1) fig ```` +We can also rotate a polygon with the same `transform` function. Here we use a +2D rotation from `Rotations.jl`, which rotates around the origin. + +````@example creating_geometry +theta = π / 4 +rotation = Rotations.Angle2d(theta); +polygon1_rotated_origin = GO.transform(p -> rotation * p, polygon1); + +fig_rotation = Figure() +ax_origin = Axis(fig_rotation[1, 1]; title = "Rotate around the origin", aspect = DataAspect()) +poly!(ax_origin, polygon1; color = (:steelblue, 0.5), strokecolor = :steelblue) +poly!(ax_origin, polygon1_rotated_origin; color = (:orange, 0.35), strokecolor = :orange) +fig_rotation +```` + +To rotate around a polygon's centroid instead, rotate each point relative to +the centroid and then shift it back. + +````@example creating_geometry +polygon1_centroid = GO.centroid(polygon1) +polygon1_rotated_centroid = GO.transform( + p -> rotation * (p .- polygon1_centroid) .+ polygon1_centroid, + polygon1, +); + +ax_centroid = Axis(fig_rotation[1, 2]; title = "Rotate around the centroid", aspect = DataAspect()) +poly!(ax_centroid, polygon1; color = (:steelblue, 0.5), strokecolor = :steelblue) +poly!(ax_centroid, polygon1_rotated_centroid; color = (:orange, 0.35), strokecolor = :orange) +fig_rotation +```` + Polygons can contain "holes". The first `LinearRing` in a polygon is the exterior, and all subsequent `LinearRing`s are treated as holes in the leading `LinearRing`. `GeoInterface` offers the `GI.getexterior(poly)` and `GI.gethole(poly)` methods to get the exterior ring and an iterable of holes, respectively. diff --git a/src/transformations/transform.jl b/src/transformations/transform.jl index 739ffad012..24b57e3ea8 100644 --- a/src/transformations/transform.jl +++ b/src/transformations/transform.jl @@ -26,9 +26,12 @@ This uses [`apply`](@ref), so will work with any geometry, vector of geometries, Apply a function `f` to all the points in `obj`. -Points will be passed to `f` as an `SVector` to allow -using CoordinateTransformations.jl and Rotations.jl -without hassle. +Points are passed to `f` as an `SVector`, so `f` can be a plain function +or a callable transform from CoordinateTransformations.jl, such as +`Translation`, `LinearMap`, or a composition of transforms. + +Because this uses [`apply`](@ref Primitive-functions) internally, it works with polygons, +multipolygons, arrays of geometries, feature collections, and tables. `SVector` is also a valid GeoInterface.jl point, so will work in all GeoInterface.jl methods. @@ -52,6 +55,37 @@ re.SVector{2, Float64}[[4.5, 3.5], [6.5, 5.5], [8.5, 7.5], [4.5, 3.5]], nothing, rraysCore.SVector{2, Float64}[[6.5, 5.5], [8.5, 7.5], [9.5, 8.5], [6.5, 5.5]], nothing, nothing)], nothing, nothing) ``` +CoordinateTransformations.jl also works directly with callable transforms like +`LinearMap`, which is handy for 2D rotation. + +```julia +julia> rotation_geom = GI.Polygon([[(0.0, 0.0), (2.0, 0.0), (2.0, 1.0), (0.0, 1.0), (0.0, 0.0)]]); + +julia> rotation = CoordinateTransformations.LinearMap([0.0 -1.0; 1.0 0.0]); + +julia> rotated = GO.transform(rotation, rotation_geom); + +julia> Tuple.(GI.getpoint(GI.getexterior(rotated))) +5-element Vector{Tuple{Float64, Float64}}: + (0.0, 0.0) + (0.0, 2.0) + (-1.0, 2.0) + (-1.0, 0.0) + (0.0, 0.0) + +julia> center = GO.centroid(rotation_geom); + +julia> rotated_centroid = GO.transform( + CoordinateTransformations.Translation(center...) ∘ + rotation ∘ + CoordinateTransformations.Translation((-).(center)...), + rotation_geom, + ); + +julia> all(GO.centroid(rotated_centroid) .≈ center) +true +``` + With Rotations.jl you need to actually multiply the Rotation by the `SVector` point, which is easy using an anonymous function. diff --git a/test/transformations/transform.jl b/test/transformations/transform.jl index 770e397866..f8ad951b66 100644 --- a/test/transformations/transform.jl +++ b/test/transformations/transform.jl @@ -6,6 +6,7 @@ using ..TestHelpers geom = GI.Polygon([GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)]), GI.LinearRing([(3, 4), (5, 6), (6, 7), (3, 4)])]) +rotation_geom = GI.Polygon([GI.LinearRing([(1.0, 1.0), (3.0, 1.0), (3.0, 2.0), (1.0, 2.0), (1.0, 1.0)])]) @testset_implementations "transform" begin translated = GI.Polygon([GI.LinearRing([[4.5, 3.5], [6.5, 5.5], [8.5, 7.5], [4.5, 3.5]]), @@ -27,3 +28,26 @@ end @test GI.is3d(geom_transformed) @test !GI.ismeasured(geom_transformed) end + +@testset_implementations "transform polygon rotation around the origin" begin + rotation = LinearMap([0.0 -1.0; 1.0 0.0]) + rotated = GO.transform(rotation, $rotation_geom) + expected_points = [ + (-1.0, 1.0), (-1.0, 3.0), (-2.0, 3.0), (-2.0, 1.0), (-1.0, 1.0), + ] + rotated_points = map(collect(GO.flatten(GI.PointTrait, rotated))) do p + (GI.x(p), GI.y(p)) + end + @test rotated_points == expected_points +end + +@testset_implementations "transform polygon rotation around centroid preserves centroid and area" begin + rotation = LinearMap([0.0 -1.0; 1.0 0.0]) + center = GO.centroid($rotation_geom) + rotated = GO.transform( + Translation(center...) ∘ rotation ∘ Translation((-).(center)...), + $rotation_geom, + ) + @test all(GO.centroid(rotated) .≈ center) + @test GO.area(rotated) ≈ GO.area($rotation_geom) +end