Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions docs/Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
32 changes: 32 additions & 0 deletions docs/src/tutorials/creating_geometry.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -113,6 +114,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.
Expand Down
40 changes: 37 additions & 3 deletions src/transformations/transform.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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) 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.
Expand All @@ -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.

Expand Down
24 changes: 24 additions & 0 deletions test/transformations/transform.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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]]),
Expand All @@ -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
Loading