From d9ff60bda7174dc3204dfabb0f8399036bd92a0b Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sun, 6 Oct 2024 14:04:32 -0700 Subject: [PATCH 01/30] Add GeometryOpsCore as an explicit dependency + fix version --- Project.toml | 2 ++ src/GeometryOps.jl | 4 ++-- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 9c067b595..a62b5f65f 100644 --- a/Project.toml +++ b/Project.toml @@ -10,6 +10,7 @@ DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df" ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110" GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f" GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" +GeometryOpsCore = "05efe853-fabf-41c8-927e-7063c8b9f013" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" SortTileRecursiveTree = "746ee33f-1797-42c2-866d-db2fce69d14d" @@ -34,6 +35,7 @@ ExactPredicates = "2.2.8" FlexiJoins = "0.1.30" GeoInterface = "1.2" GeometryBasics = "0.4.7" +GeometryOpsCore = "=0.1.1" LibGEOS = "0.9.2" LinearAlgebra = "1" Proj = "1" diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index f524fd87d..92f94ef1f 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -2,8 +2,8 @@ module GeometryOps -include("../GeometryOpsCore/src/GeometryOpsCore.jl") # TODO: replace this with `using GeometryOpsCore` -import .GeometryOpsCore +import GeometryOpsCore +for name in Base.setdiff( for name in setdiff(names(GeometryOpsCore, all = true), (:eval, :var"#eval", :include, :var"#include")) # Import all symbols from GeometryOpsCore @eval import .GeometryOpsCore: $name From f50c82a0a2137fa8be6534afe550f1b2fb949f6e Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sun, 6 Oct 2024 14:04:42 -0700 Subject: [PATCH 02/30] Fix imports to work on nightly --- src/GeometryOps.jl | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index 92f94ef1f..8373d470e 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -4,9 +4,14 @@ module GeometryOps import GeometryOpsCore for name in Base.setdiff( -for name in setdiff(names(GeometryOpsCore, all = true), (:eval, :var"#eval", :include, :var"#include")) + Base.union( + names(GeometryOpsCore), + (:flatten, :reconstruct, :rebuild, :unwrap, :APPLY_KEYWORDS, :THREADED_KEYWORD, :CRS_KEYWORD, :CALC_EXTENT_KEYWORD,) + ), + (:eval, :include, :var"#eval", :var"#include"), + ) # Import all symbols from GeometryOpsCore - @eval import .GeometryOpsCore: $name + @eval import GeometryOpsCore: $name # Re-export all exported symbols if Base.isexported(GeometryOpsCore, name) @eval export $name From 324fccbaf208e85863e22d6c0d1d8e5a1f2863e0 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sun, 6 Oct 2024 14:06:23 -0700 Subject: [PATCH 03/30] Add nightly job --- .github/workflows/CI.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 84a14a375..c293da4d7 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -20,7 +20,7 @@ jobs: version: - '1.9' - '1' - # - 'nightly' + - 'nightly' os: - ubuntu-latest arch: From 36bb92ca871e13bb94abc6ea2f6420fecb45d6af Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sun, 6 Oct 2024 14:06:34 -0700 Subject: [PATCH 04/30] Beautify --- src/GeometryOps.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index 8373d470e..b5c8ed486 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -6,7 +6,10 @@ import GeometryOpsCore for name in Base.setdiff( Base.union( names(GeometryOpsCore), - (:flatten, :reconstruct, :rebuild, :unwrap, :APPLY_KEYWORDS, :THREADED_KEYWORD, :CRS_KEYWORD, :CALC_EXTENT_KEYWORD,) + ( + :flatten, :reconstruct, :rebuild, :unwrap, + :APPLY_KEYWORDS, :THREADED_KEYWORD, :CRS_KEYWORD, :CALC_EXTENT_KEYWORD, + ) ), (:eval, :include, :var"#eval", :var"#include"), ) From e40db7d27071fc2efdc2e0389180b414e8f37b40 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sun, 6 Oct 2024 14:11:37 -0700 Subject: [PATCH 05/30] import linearring as well --- src/GeometryOps.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index b5c8ed486..c4d16ff4c 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -7,7 +7,7 @@ for name in Base.setdiff( Base.union( names(GeometryOpsCore), ( - :flatten, :reconstruct, :rebuild, :unwrap, + :flatten, :reconstruct, :rebuild, :unwrap, :_linearring, :APPLY_KEYWORDS, :THREADED_KEYWORD, :CRS_KEYWORD, :CALC_EXTENT_KEYWORD, ) ), From 200946321a0fd09d8060cf7e3993f50c21f537d9 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sun, 6 Oct 2024 22:02:37 -0700 Subject: [PATCH 06/30] Simplify and make imports explicit --- src/GeometryOps.jl | 25 ++++++++----------------- 1 file changed, 8 insertions(+), 17 deletions(-) diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index c4d16ff4c..90f0da2e1 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -3,23 +3,14 @@ module GeometryOps import GeometryOpsCore -for name in Base.setdiff( - Base.union( - names(GeometryOpsCore), - ( - :flatten, :reconstruct, :rebuild, :unwrap, :_linearring, - :APPLY_KEYWORDS, :THREADED_KEYWORD, :CRS_KEYWORD, :CALC_EXTENT_KEYWORD, - ) - ), - (:eval, :include, :var"#eval", :var"#include"), - ) - # Import all symbols from GeometryOpsCore - @eval import GeometryOpsCore: $name - # Re-export all exported symbols - if Base.isexported(GeometryOpsCore, name) - @eval export $name - end -end +import GeometryOpsCore: + Manifold, Planar, Spherical, Geodesic, + BoolsAsTypes, _True, _False, _booltype, + apply, applyreduce, + flatten, reconstruct, rebuild, unwrap, _linearring, + APPLY_KEYWORDS, THREADED_KEYWORD, CRS_KEYWORD, CALC_EXTENT_KEYWORD + +export Manifold, Planar, Spherical, Geodesic, apply, applyreduce, flatten, reconstruct, rebuild, unwrap using GeoInterface using GeometryBasics From 611a726b82e4d1dfd3b9f24c36b651d5c81cb584 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sun, 6 Oct 2024 22:28:13 -0700 Subject: [PATCH 07/30] impex TraitTarget + make sure GO and GOC are both local --- .github/workflows/CI.yml | 4 +++- src/GeometryOps.jl | 3 ++- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index c293da4d7..4bf45fd31 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -33,6 +33,8 @@ jobs: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - uses: julia-actions/cache@v1 + - name: Dev GeometryOpsCore + run: julia --project=. -e 'using Pkg; Pkg.dev(joinpath(".", "GeometryOpsCore"))' - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 @@ -58,7 +60,7 @@ jobs: with: version: '1' - name: Build and add versions - run: julia --project=docs -e 'using Pkg; Pkg.add([PackageSpec(path = "."), PackageSpec(name = "GeoMakie", rev = "master")])' + run: julia --project=docs -e 'using Pkg; Pkg.add([PackageSpec(path = "."), PackageSpec(path = joinpath(".", "GeometryOpsCore")), PackageSpec(name = "GeoMakie", rev = "master")])' - uses: julia-actions/julia-docdeploy@v1 with: install-package: false diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index 90f0da2e1..e0cd484bd 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -4,13 +4,14 @@ module GeometryOps import GeometryOpsCore import GeometryOpsCore: + TraitTarget, Manifold, Planar, Spherical, Geodesic, BoolsAsTypes, _True, _False, _booltype, apply, applyreduce, flatten, reconstruct, rebuild, unwrap, _linearring, APPLY_KEYWORDS, THREADED_KEYWORD, CRS_KEYWORD, CALC_EXTENT_KEYWORD -export Manifold, Planar, Spherical, Geodesic, apply, applyreduce, flatten, reconstruct, rebuild, unwrap +export TraitTarget, Manifold, Planar, Spherical, Geodesic, apply, applyreduce, flatten, reconstruct, rebuild, unwrap using GeoInterface using GeometryBasics From 86f3b8766fa6bd8cec23c47aabb99c82e6ea41e7 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sun, 6 Oct 2024 22:35:39 -0700 Subject: [PATCH 08/30] Use add not dev in CI --- .github/workflows/CI.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 4bf45fd31..8e3271b9b 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -33,8 +33,8 @@ jobs: version: ${{ matrix.version }} arch: ${{ matrix.arch }} - uses: julia-actions/cache@v1 - - name: Dev GeometryOpsCore - run: julia --project=. -e 'using Pkg; Pkg.dev(joinpath(".", "GeometryOpsCore"))' + - name: Dev GeometryOpsCore` + run: julia --project=. -e 'using Pkg; Pkg.add(PackageSpec(; path = joinpath(".", "GeometryOpsCore")))' - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 From 73f486b12974dbc48da66f7070c7bf6bd0a743fc Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sun, 6 Oct 2024 22:36:01 -0700 Subject: [PATCH 09/30] Adapt segmentize to the new way of doing things --- ext/GeometryOpsProjExt/segmentize.jl | 25 +++++++++++++++--- src/transformations/segmentize.jl | 39 ++++++++++++++++++---------- test/transformations/segmentize.jl | 21 +++++++-------- 3 files changed, 57 insertions(+), 28 deletions(-) diff --git a/ext/GeometryOpsProjExt/segmentize.jl b/ext/GeometryOpsProjExt/segmentize.jl index 308c97188..40935ec44 100644 --- a/ext/GeometryOpsProjExt/segmentize.jl +++ b/ext/GeometryOpsProjExt/segmentize.jl @@ -7,14 +7,31 @@ function GeometryOps.GeodesicSegments(; max_distance, equatorial_radius::Real=63 return GeometryOps.GeodesicSegments{Proj.geod_geodesic}(geodesic, max_distance) end +# This is the same method as in `transformations/segmentize.jl`, +# but it constructs a Proj geodesic line every time. +# Maybe this should be better... +function _segmentize(method::Geodesic, geom, T::Union{GI.LineStringTrait, GI.LinearRingTrait}; max_distance) + proj_geodesic = Proj.geod_geodesic(method.equatorial_radius, method.flattening) + first_coord = GI.getpoint(geom, 1) + x1, y1 = GI.x(first_coord), GI.y(first_coord) + new_coords = NTuple{2, Float64}[] + sizehint!(new_coords, GI.npoint(geom)) + push!(new_coords, (x1, y1)) + for coord in Iterators.drop(GI.getpoint(geom), 1) + x2, y2 = GI.x(coord), GI.y(coord) + _fill_linear_kernel!(method, new_coords, x1, y1, x2, y2; max_distance, proj_geodesic) + x1, y1 = x2, y2 + end + return rebuild(geom, new_coords) +end -function GeometryOps._fill_linear_kernel!(method::GeodesicSegments{Proj.geod_geodesic}, new_coords::Vector, x1, y1, x2, y2) - geod_line = Proj.geod_inverseline(method.geodesic, y1, x1, y2, x2) +function GeometryOps._fill_linear_kernel!(method::Geodesic, new_coords::Vector, x1, y1, x2, y2; max_distance, proj_geodesic) + geod_line = Proj.geod_inverseline(proj_geodesic, y1, x1, y2, x2) # This is the distance in meters computed between the two points. # It's `s13` because `geod_inverseline` sets point 3 to the second input point. distance = geod_line.s13 - if distance > method.max_distance - n_segments = ceil(Int, distance / method.max_distance) + if distance > max_distance + n_segments = ceil(Int, distance / max_distance) for i in 1:(n_segments - 1) y, x, _ = Proj.geod_position(geod_line, i / n_segments * distance) push!(new_coords, (x, y)) diff --git a/src/transformations/segmentize.jl b/src/transformations/segmentize.jl index 623637bc5..ef653334b 100644 --- a/src/transformations/segmentize.jl +++ b/src/transformations/segmentize.jl @@ -146,7 +146,7 @@ end # Add an error hint for GeodesicSegments if Proj is not loaded! function _geodesic_segments_error_hinter(io, exc, argtypes, kwargs) if isnothing(Base.get_extension(GeometryOps, :GeometryOpsProjExt)) && exc.f == GeodesicSegments - print(io, "\n\nThe `GeodesicSegments` method requires the Proj.jl package to be explicitly loaded.\n") + print(io, "\n\nThe `Geodesic` method requires the Proj.jl package to be explicitly loaded.\n") print(io, "You can do this by simply typing ") printstyled(io, "using Proj"; color = :cyan, bold = true) println(io, " in your REPL, \nor otherwise loading Proj.jl via using or import.") @@ -156,33 +156,46 @@ end # ## Implementation """ - segmentize([method = LinearSegments()], geom; max_distance::Real, threaded) + segmentize([method = Planar()], geom; max_distance::Real, threaded) Segmentize a geometry by adding extra vertices to the geometry so that no segment is longer than a given distance. This is useful for plotting geometries with a limited number of vertices, or for ensuring that a geometry is not too "coarse" for a given application. ## Arguments -- `method::SegmentizeMethod = LinearSegments()`: The method to use for segmentizing the geometry. At the moment, only [`LinearSegments`](@ref) and [`GeodesicSegments`](@ref) are available. -- `geom`: The geometry to segmentize. Must be a `LineString`, `LinearRing`, or greater in complexity. -- `max_distance::Real`: The maximum distance, **in the input space**, between vertices in the geometry. Only used if you don't explicitly pass a `method`. +- `method::Manifold = Planar()`: The method to use for segmentizing the geometry. At the moment, only [`Planar`](@ref) (assumes a flat plane) and [`Geodesic`](@ref) (assumes geometry on the ellipsoidal Earth and uses Vincenty's formulae) are available. +- `geom`: The geometry to segmentize. Must be a `LineString`, `LinearRing`, `Polygon`, `MultiPolygon`, or `GeometryCollection`, or some vector or table of those. +- `max_distance::Real`: The maximum distance between vertices in the geometry. **Beware: for `Planar`, this is in the units of the geometry, but for `Geodesic` and `Spherical` it's in units of the radius of the sphere.** Returns a geometry of similar type to the input geometry, but resampled. """ function segmentize(geom; max_distance, threaded::Union{Bool, BoolsAsTypes} = _False()) - return segmentize(LinearSegments(; max_distance), geom; threaded = _booltype(threaded)) + return segmentize(Linear(), geom; threaded = _booltype(threaded)) end -function segmentize(method::SegmentizeMethod, geom; threaded::Union{Bool, BoolsAsTypes} = _False()) - @assert method.max_distance > 0 "`max_distance` should be positive and nonzero! Found $(method.max_distance)." + +# allow three-arg method as well, just in case +segmentize(geom, max_distance::Real; threaded = _False()) = segmentize(Linear(), geom, max_distance; threaded) +segmentize(method::Manifold, geom, max_distance::Real; threaded = _False()) = segmentize(Linear(), geom; max_distance, threaded) + +# generic implementation +function segmentize(method::Manifold, geom; max_distance, threaded::Union{Bool, BoolsAsTypes} = _False()) + @assert max_distance > 0 "`max_distance` should be positive and nonzero! Found $(method.max_distance)." segmentize_function = Base.Fix1(_segmentize, method) return apply(segmentize_function, TraitTarget(GI.LinearRingTrait(), GI.LineStringTrait()), geom; threaded) end +function segmentize(method::SegmentizeMethod, geom; threaded::Union{Bool, BoolsAsTypes} = _False()) + @warn "`segmentize(method::$(typeof(method)), geom) is deprecated; use `segmentize($(method isa LinearSegments ? "Linear()" : "Geodesic()"), geom; max_distance, threaded) instead!" maxlog=3 + @assert method.max_distance > 0 "`max_distance` should be positive and nonzero! Found $(method.max_distance)." + new_method = method isa LinearSegments ? Linear() : Geodesic() + segmentize(new_method, geom; max_distance = method.max_distance, threaded) +end + _segmentize(method, geom) = _segmentize(method, geom, GI.trait(geom)) #= This is a method which performs the common functionality for both linear and geodesic algorithms, and calls out to the "kernel" function which we've defined per linesegment. =# -function _segmentize(method::Union{LinearSegments, GeodesicSegments}, geom, T::Union{GI.LineStringTrait, GI.LinearRingTrait}) +function _segmentize(method::Union{Linear, Spherical}, geom, T::Union{GI.LineStringTrait, GI.LinearRingTrait}; max_distance) first_coord = GI.getpoint(geom, 1) x1, y1 = GI.x(first_coord), GI.y(first_coord) new_coords = NTuple{2, Float64}[] @@ -190,17 +203,17 @@ function _segmentize(method::Union{LinearSegments, GeodesicSegments}, geom, T::U push!(new_coords, (x1, y1)) for coord in Iterators.drop(GI.getpoint(geom), 1) x2, y2 = GI.x(coord), GI.y(coord) - _fill_linear_kernel!(method, new_coords, x1, y1, x2, y2) + _fill_linear_kernel!(method, new_coords, x1, y1, x2, y2; max_distance) x1, y1 = x2, y2 end return rebuild(geom, new_coords) end -function _fill_linear_kernel!(method::LinearSegments, new_coords::Vector, x1, y1, x2, y2) +function _fill_linear_kernel!(::Linear, new_coords::Vector, x1, y1, x2, y2; max_distance) dx, dy = x2 - x1, y2 - y1 distance = hypot(dx, dy) # this is a more stable way to compute the Euclidean distance - if distance > method.max_distance - n_segments = ceil(Int, distance / method.max_distance) + if distance > max_distance + n_segments = ceil(Int, distance / max_distance) for i in 1:(n_segments - 1) t = i / n_segments push!(new_coords, (x1 + t * dx, y1 + t * dy)) diff --git a/test/transformations/segmentize.jl b/test/transformations/segmentize.jl index 5095c4c84..9746a5d0a 100644 --- a/test/transformations/segmentize.jl +++ b/test/transformations/segmentize.jl @@ -31,21 +31,20 @@ using ..TestHelpers end @testset_implementations "GeodesicSegments" begin - @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), $ls) isa GI.LineString + @test GO.segmentize(GO.Geodesic(), $ls; max_distance = 0.5*900) isa GI.LineString if GI.trait($lr) isa GI.LinearRingTrait - @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), $lr) isa GI.LinearRing + @test GO.segmentize(GO.Geodesic(), $lr; max_distance = 0.5*900) isa GI.LinearRing end # Test that linear rings are closed after segmentization - segmentized_ring = GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), $lr) + segmentized_ring = GO.segmentize(GO.Geodesic(), $lr; max_distance = 0.5*900) @test GI.getpoint(segmentized_ring, 1) == GI.getpoint(segmentized_ring, GI.npoint(segmentized_ring)) - - @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), $p) isa GI.Polygon - @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), $mp) isa GI.MultiPolygon - @test GI.ngeom(GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), $mp)) == 3 + @test GO.segmentize(GO.Geodesic(), $p; max_distance = 0.5*900) isa GI.Polygon + @test GO.segmentize(GO.Geodesic(), $mp; max_distance = 0.5*900) isa GI.MultiPolygon + @test GI.ngeom(GO.segmentize(GO.Geodesic(), $mp; max_distance = 0.5*900)) == 3 # Now test multilinestrings - @test GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), $mls) isa GI.MultiLineString - @test GI.ngeom(GO.segmentize(GO.GeodesicSegments(; max_distance = 0.5*900), $mls)) == 3 + @test GO.segmentize(GO.Geodesic(), $mls; max_distance = 0.5*900) isa GI.MultiLineString + @test GI.ngeom(GO.segmentize(GO.Geodesic(), $mls; max_distance = 0.5*900)) == 3 end end @@ -55,7 +54,7 @@ lr = GI.LinearRing([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)]) ct = GO.centroid($lr) ar = GO.area($lr) for max_distance in exp10.(LinRange(log10(0.01), log10(1), 10)) - segmentized = GO.segmentize(GO.LinearSegments(; max_distance), $lr) + segmentized = GO.segmentize(GO.Linear(), $lr; max_distance) @test all(GO.centroid(segmentized) .≈ ct) @test GO.area(segmentized) ≈ ar end @@ -64,6 +63,6 @@ end lr = GI.LinearRing([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)]) @testset_implementations "GeodesicSegments" begin for max_distance in exp10.(LinRange(log10(0.01), log10(1), 10)) .* 900 - @test_nowarn segmentized = GO.segmentize(GO.GeodesicSegments(; max_distance), $lr) + @test_nowarn segmentized = GO.segmentize(GO.Geodesic(), $lr; max_distance) end end From 1f5613a3e4eef4baea4d230e71849f4d7ab7160f Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sun, 6 Oct 2024 22:42:29 -0700 Subject: [PATCH 10/30] try and fix ci --- .github/workflows/CI.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 8e3271b9b..fab30bff7 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -34,7 +34,7 @@ jobs: arch: ${{ matrix.arch }} - uses: julia-actions/cache@v1 - name: Dev GeometryOpsCore` - run: julia --project=. -e 'using Pkg; Pkg.add(PackageSpec(; path = joinpath(".", "GeometryOpsCore")))' + run: julia --project=. -e 'using Pkg; Pkg.develop(; path = joinpath(".", "GeometryOpsCore"))' - uses: julia-actions/julia-buildpkg@v1 - uses: julia-actions/julia-runtest@v1 - uses: julia-actions/julia-processcoverage@v1 From b07808fe2de0514f03d22086f8e9b0b18bf7e9bc Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sun, 6 Oct 2024 22:43:19 -0700 Subject: [PATCH 11/30] devadd in docs too --- .github/workflows/CI.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index fab30bff7..cbbb3e165 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -60,7 +60,7 @@ jobs: with: version: '1' - name: Build and add versions - run: julia --project=docs -e 'using Pkg; Pkg.add([PackageSpec(path = "."), PackageSpec(path = joinpath(".", "GeometryOpsCore")), PackageSpec(name = "GeoMakie", rev = "master")])' + run: julia --project=docs -e 'using Pkg; Pkg.develop([PackageSpec(path = "."), PackageSpec(path = joinpath(".", "GeometryOpsCore"))]); Pkg.add([PackageSpec(name = "GeoMakie", rev = "master")])' - uses: julia-actions/julia-docdeploy@v1 with: install-package: false From 7ea0676ae098ac479e53f1528cdaf69f91606555 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sun, 6 Oct 2024 22:49:04 -0700 Subject: [PATCH 12/30] linear -> planar --- src/transformations/segmentize.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/transformations/segmentize.jl b/src/transformations/segmentize.jl index ef653334b..d4d147e87 100644 --- a/src/transformations/segmentize.jl +++ b/src/transformations/segmentize.jl @@ -169,12 +169,12 @@ This is useful for plotting geometries with a limited number of vertices, or for Returns a geometry of similar type to the input geometry, but resampled. """ function segmentize(geom; max_distance, threaded::Union{Bool, BoolsAsTypes} = _False()) - return segmentize(Linear(), geom; threaded = _booltype(threaded)) + return segmentize(Planar(), geom; threaded = _booltype(threaded)) end # allow three-arg method as well, just in case -segmentize(geom, max_distance::Real; threaded = _False()) = segmentize(Linear(), geom, max_distance; threaded) -segmentize(method::Manifold, geom, max_distance::Real; threaded = _False()) = segmentize(Linear(), geom; max_distance, threaded) +segmentize(geom, max_distance::Real; threaded = _False()) = segmentize(Planar(), geom, max_distance; threaded) +segmentize(method::Manifold, geom, max_distance::Real; threaded = _False()) = segmentize(Planar(), geom; max_distance, threaded) # generic implementation function segmentize(method::Manifold, geom; max_distance, threaded::Union{Bool, BoolsAsTypes} = _False()) @@ -184,9 +184,9 @@ function segmentize(method::Manifold, geom; max_distance, threaded::Union{Bool, end function segmentize(method::SegmentizeMethod, geom; threaded::Union{Bool, BoolsAsTypes} = _False()) - @warn "`segmentize(method::$(typeof(method)), geom) is deprecated; use `segmentize($(method isa LinearSegments ? "Linear()" : "Geodesic()"), geom; max_distance, threaded) instead!" maxlog=3 + @warn "`segmentize(method::$(typeof(method)), geom) is deprecated; use `segmentize($(method isa LinearSegments ? "Planar()" : "Geodesic()"), geom; max_distance, threaded) instead!" maxlog=3 @assert method.max_distance > 0 "`max_distance` should be positive and nonzero! Found $(method.max_distance)." - new_method = method isa LinearSegments ? Linear() : Geodesic() + new_method = method isa LinearSegments ? Planar() : Geodesic() segmentize(new_method, geom; max_distance = method.max_distance, threaded) end From 50cdbef6665aa82501d402aa5b0e66482703cf22 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sun, 6 Oct 2024 23:04:32 -0700 Subject: [PATCH 13/30] more linear -> planar --- src/transformations/segmentize.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/transformations/segmentize.jl b/src/transformations/segmentize.jl index d4d147e87..a14f6ac7a 100644 --- a/src/transformations/segmentize.jl +++ b/src/transformations/segmentize.jl @@ -195,7 +195,7 @@ _segmentize(method, geom) = _segmentize(method, geom, GI.trait(geom)) This is a method which performs the common functionality for both linear and geodesic algorithms, and calls out to the "kernel" function which we've defined per linesegment. =# -function _segmentize(method::Union{Linear, Spherical}, geom, T::Union{GI.LineStringTrait, GI.LinearRingTrait}; max_distance) +function _segmentize(method::Union{Planar, Spherical}, geom, T::Union{GI.LineStringTrait, GI.LinearRingTrait}; max_distance) first_coord = GI.getpoint(geom, 1) x1, y1 = GI.x(first_coord), GI.y(first_coord) new_coords = NTuple{2, Float64}[] @@ -209,7 +209,7 @@ function _segmentize(method::Union{Linear, Spherical}, geom, T::Union{GI.LineStr return rebuild(geom, new_coords) end -function _fill_linear_kernel!(::Linear, new_coords::Vector, x1, y1, x2, y2; max_distance) +function _fill_linear_kernel!(::Planar, new_coords::Vector, x1, y1, x2, y2; max_distance) dx, dy = x2 - x1, y2 - y1 distance = hypot(dx, dy) # this is a more stable way to compute the Euclidean distance if distance > max_distance From 91ee1cad2e896e30e16f41e6f99765275b044786 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 7 Oct 2024 20:17:16 -0700 Subject: [PATCH 14/30] Core: fix typos, bump version --- GeometryOpsCore/Project.toml | 2 +- GeometryOpsCore/src/types.jl | 4 +++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/GeometryOpsCore/Project.toml b/GeometryOpsCore/Project.toml index b364aa25a..e3bc961d8 100644 --- a/GeometryOpsCore/Project.toml +++ b/GeometryOpsCore/Project.toml @@ -1,7 +1,7 @@ name = "GeometryOpsCore" uuid = "05efe853-fabf-41c8-927e-7063c8b9f013" authors = ["Anshul Singhvi ", "Rafael Schouten ", "Skylar Gering ", "and contributors"] -version = "0.1.1" +version = "0.1.2" [deps] DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" diff --git a/GeometryOpsCore/src/types.jl b/GeometryOpsCore/src/types.jl index 5bc962a71..fbe44fcef 100644 --- a/GeometryOpsCore/src/types.jl +++ b/GeometryOpsCore/src/types.jl @@ -51,6 +51,8 @@ end A spherical manifold means that the geometry is on the 3-sphere (but is represented by 2-D longitude and latitude). +By default, the radius is defined as the mean radius of the Earth, `6371008.8 m`. + ## Extended help !!! note @@ -74,7 +76,7 @@ and `inv_flattening` (``1/f``). Usually, this is only relevant for area and segmentization calculations. It becomes more relevant as one grows closer to the poles (or equator). """ Base.@kwdef struct Geodesic{T} <: Manifold - semimajor_axis::T = 6378137,0 + semimajor_axis::T = 6378137.0 inv_flattening::T = 298.257223563 end From 739c55b3bd5db3bcae5965c558818d3347494e7a Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 7 Oct 2024 20:17:20 -0700 Subject: [PATCH 15/30] Bump Core compat --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index a62b5f65f..39fbe34e3 100644 --- a/Project.toml +++ b/Project.toml @@ -35,7 +35,7 @@ ExactPredicates = "2.2.8" FlexiJoins = "0.1.30" GeoInterface = "1.2" GeometryBasics = "0.4.7" -GeometryOpsCore = "=0.1.1" +GeometryOpsCore = "=0.1.2" LibGEOS = "0.9.2" LinearAlgebra = "1" Proj = "1" From a41c25714a60209cf650827420d8b508fb6f4686 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 7 Oct 2024 20:31:26 -0700 Subject: [PATCH 16/30] fix default segmentize not passing on max_distance --- src/transformations/segmentize.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/transformations/segmentize.jl b/src/transformations/segmentize.jl index a14f6ac7a..fdd8d1985 100644 --- a/src/transformations/segmentize.jl +++ b/src/transformations/segmentize.jl @@ -169,7 +169,7 @@ This is useful for plotting geometries with a limited number of vertices, or for Returns a geometry of similar type to the input geometry, but resampled. """ function segmentize(geom; max_distance, threaded::Union{Bool, BoolsAsTypes} = _False()) - return segmentize(Planar(), geom; threaded = _booltype(threaded)) + return segmentize(Planar(), geom; max_distance, threaded = _booltype(threaded)) end # allow three-arg method as well, just in case From 53936f1af7f695aa21e1ff2e05c31a3f7ae56ddf Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 7 Oct 2024 20:43:39 -0700 Subject: [PATCH 17/30] fix geodesic again --- ext/GeometryOpsProjExt/segmentize.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ext/GeometryOpsProjExt/segmentize.jl b/ext/GeometryOpsProjExt/segmentize.jl index 40935ec44..e933f17a2 100644 --- a/ext/GeometryOpsProjExt/segmentize.jl +++ b/ext/GeometryOpsProjExt/segmentize.jl @@ -1,6 +1,6 @@ # This holds the `segmentize` geodesic functionality. -import GeometryOps: GeodesicSegments, _fill_linear_kernel! +import GeometryOps: GeodesicSegments, _segmentize, _fill_linear_kernel! import Proj function GeometryOps.GeodesicSegments(; max_distance, equatorial_radius::Real=6378137, flattening::Real=1/298.257223563, geodesic::Proj.geod_geodesic = Proj.geod_geodesic(equatorial_radius, flattening)) @@ -10,7 +10,7 @@ end # This is the same method as in `transformations/segmentize.jl`, # but it constructs a Proj geodesic line every time. # Maybe this should be better... -function _segmentize(method::Geodesic, geom, T::Union{GI.LineStringTrait, GI.LinearRingTrait}; max_distance) +function _segmentize(method::Geodesic, geom, ::Union{GI.LineStringTrait, GI.LinearRingTrait}; max_distance) proj_geodesic = Proj.geod_geodesic(method.equatorial_radius, method.flattening) first_coord = GI.getpoint(geom, 1) x1, y1 = GI.x(first_coord), GI.y(first_coord) From 60593d6fc3dbd8084d8a90e5a2e9e000a250e40a Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 7 Oct 2024 20:50:57 -0700 Subject: [PATCH 18/30] fix all segmentizes --- src/transformations/segmentize.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/transformations/segmentize.jl b/src/transformations/segmentize.jl index fdd8d1985..32f42d2ff 100644 --- a/src/transformations/segmentize.jl +++ b/src/transformations/segmentize.jl @@ -179,7 +179,7 @@ segmentize(method::Manifold, geom, max_distance::Real; threaded = _False()) = se # generic implementation function segmentize(method::Manifold, geom; max_distance, threaded::Union{Bool, BoolsAsTypes} = _False()) @assert max_distance > 0 "`max_distance` should be positive and nonzero! Found $(method.max_distance)." - segmentize_function = Base.Fix1(_segmentize, method) + segmentize_function(x) = _segmentize(method, x; max_distance) return apply(segmentize_function, TraitTarget(GI.LinearRingTrait(), GI.LineStringTrait()), geom; threaded) end From 3841ec8e9d184aa7e28d00bb2a52297d5696a4fa Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 7 Oct 2024 21:12:29 -0700 Subject: [PATCH 19/30] +fix --- src/transformations/segmentize.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/transformations/segmentize.jl b/src/transformations/segmentize.jl index 32f42d2ff..df5882406 100644 --- a/src/transformations/segmentize.jl +++ b/src/transformations/segmentize.jl @@ -179,8 +179,8 @@ segmentize(method::Manifold, geom, max_distance::Real; threaded = _False()) = se # generic implementation function segmentize(method::Manifold, geom; max_distance, threaded::Union{Bool, BoolsAsTypes} = _False()) @assert max_distance > 0 "`max_distance` should be positive and nonzero! Found $(method.max_distance)." - segmentize_function(x) = _segmentize(method, x; max_distance) - return apply(segmentize_function, TraitTarget(GI.LinearRingTrait(), GI.LineStringTrait()), geom; threaded) + _segmentize_function(geom) = _segmentize(method, geom, GI.trait(geom); max_distance) + return apply(_segmentize_function, TraitTarget(GI.LinearRingTrait(), GI.LineStringTrait()), geom; threaded) end function segmentize(method::SegmentizeMethod, geom; threaded::Union{Bool, BoolsAsTypes} = _False()) From 617764061f1d882505e3fdc36d06848cae217c2e Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sun, 10 Nov 2024 23:03:30 -0500 Subject: [PATCH 20/30] Fix method error with geodesic --- ext/GeometryOpsProjExt/segmentize.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/GeometryOpsProjExt/segmentize.jl b/ext/GeometryOpsProjExt/segmentize.jl index e933f17a2..a7783ea02 100644 --- a/ext/GeometryOpsProjExt/segmentize.jl +++ b/ext/GeometryOpsProjExt/segmentize.jl @@ -11,7 +11,7 @@ end # but it constructs a Proj geodesic line every time. # Maybe this should be better... function _segmentize(method::Geodesic, geom, ::Union{GI.LineStringTrait, GI.LinearRingTrait}; max_distance) - proj_geodesic = Proj.geod_geodesic(method.equatorial_radius, method.flattening) + proj_geodesic = Proj.geod_geodesic(method.semimajor_axis #= same thing as equatorial radius =#, 1/method.inv_flattening) first_coord = GI.getpoint(geom, 1) x1, y1 = GI.x(first_coord), GI.y(first_coord) new_coords = NTuple{2, Float64}[] From a38f74aa9860a78a202c035bb54139d2cd1cb7e4 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 11 Nov 2024 17:10:21 -0500 Subject: [PATCH 21/30] change Linear to Planar --- test/transformations/segmentize.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/transformations/segmentize.jl b/test/transformations/segmentize.jl index 9746a5d0a..10d07c7f3 100644 --- a/test/transformations/segmentize.jl +++ b/test/transformations/segmentize.jl @@ -50,18 +50,18 @@ using ..TestHelpers end lr = GI.LinearRing([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)]) -@testset_implementations "LinearSegments" begin +@testset_implementations "Planar" begin ct = GO.centroid($lr) ar = GO.area($lr) for max_distance in exp10.(LinRange(log10(0.01), log10(1), 10)) - segmentized = GO.segmentize(GO.Linear(), $lr; max_distance) + segmentized = GO.segmentize(GO.Planar(), $lr; max_distance) @test all(GO.centroid(segmentized) .≈ ct) @test GO.area(segmentized) ≈ ar end end lr = GI.LinearRing([(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)]) -@testset_implementations "GeodesicSegments" begin +@testset_implementations "Geodesic" begin for max_distance in exp10.(LinRange(log10(0.01), log10(1), 10)) .* 900 @test_nowarn segmentized = GO.segmentize(GO.Geodesic(), $lr; max_distance) end From 2c7902871c016f0036bfa51c9308348366d41a1f Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Tue, 12 Nov 2024 13:51:28 -0500 Subject: [PATCH 22/30] load booltypes --- ext/GeometryOpsLibGEOSExt/GeometryOpsLibGEOSExt.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ext/GeometryOpsLibGEOSExt/GeometryOpsLibGEOSExt.jl b/ext/GeometryOpsLibGEOSExt/GeometryOpsLibGEOSExt.jl index c8cbfc9f8..1d0845ed0 100644 --- a/ext/GeometryOpsLibGEOSExt/GeometryOpsLibGEOSExt.jl +++ b/ext/GeometryOpsLibGEOSExt/GeometryOpsLibGEOSExt.jl @@ -3,7 +3,7 @@ module GeometryOpsLibGEOSExt import GeometryOps as GO, LibGEOS as LG import GeoInterface as GI -import GeometryOps: GEOS, enforce +import GeometryOps: GEOS, enforce, _True, _False, _booltype using GeometryOps # The filter statement is required because in Julia, each module has its own versions of these From 39b51f9d7e8db06c99feea5846aba8f9cd970be3 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 25 Nov 2024 17:15:08 -0500 Subject: [PATCH 23/30] Remove useless functions from not_implemented_yet (#237) --- src/not_implemented_yet.jl | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/not_implemented_yet.jl b/src/not_implemented_yet.jl index 45e23ef55..4e8cb8208 100644 --- a/src/not_implemented_yet.jl +++ b/src/not_implemented_yet.jl @@ -3,7 +3,4 @@ # Some of them may have implementations in LibGEOS which we can use # via an extension, but there is no native-Julia implementation for them. -function symdifference end -function buffer end -function convexhull end -function concavehull end \ No newline at end of file +function symdifference end \ No newline at end of file From ae6c92c25a301f1bdea6a2532f8f03865b962f1a Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 9 Jan 2025 10:23:53 -0100 Subject: [PATCH 24/30] Fix the inference error by removing assume_effects (#245) --- GeometryOpsCore/src/apply.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/GeometryOpsCore/src/apply.jl b/GeometryOpsCore/src/apply.jl index d7713195d..06e2ae545 100644 --- a/GeometryOpsCore/src/apply.jl +++ b/GeometryOpsCore/src/apply.jl @@ -343,9 +343,12 @@ using Base.Threads: nthreads, @threads, @spawn return mapreduce(fetch, vcat, tasks) end #= -Here we use the compiler directive `@assume_effects :foldable` to force the compiler +Here we used to use the compiler directive `@assume_effects :foldable` to force the compiler to lookup through the closure. This alone makes e.g. `flip` 2.5x faster! + +But it caused inference to fail, so we've removed it. No effect on runtime so far as we can tell, +at least in Julia 1.11. =# -Base.@assume_effects :foldable @inline function _maptasks(f::F, taskrange, threaded::_False)::Vector where F +@inline function _maptasks(f::F, taskrange, threaded::_False)::Vector where F map(f, taskrange) end \ No newline at end of file From f526f54431dbc2d001d9ca645ac3c959254a4fa1 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Fri, 10 Jan 2025 11:49:17 -0100 Subject: [PATCH 25/30] Check the dimensions of child geoms when rebuilding + add forcexy, forcexyz (#239) --- GeometryOpsCore/src/other_primitives.jl | 12 ++++----- src/GeometryOps.jl | 1 + src/transformations/forcedims.jl | 36 +++++++++++++++++++++++++ test/runtests.jl | 1 + test/transformations/forcedims.jl | 35 ++++++++++++++++++++++++ test/transformations/transform.jl | 14 ++++++++++ 6 files changed, 93 insertions(+), 6 deletions(-) create mode 100644 src/transformations/forcedims.jl create mode 100644 test/transformations/forcedims.jl diff --git a/GeometryOpsCore/src/other_primitives.jl b/GeometryOpsCore/src/other_primitives.jl index b0ac8c6b7..c165f525e 100644 --- a/GeometryOpsCore/src/other_primitives.jl +++ b/GeometryOpsCore/src/other_primitives.jl @@ -159,10 +159,10 @@ on geometries from other packages and specify how to rebuild them. rebuild(geom, child_geoms; kw...) = rebuild(GI.trait(geom), geom, child_geoms; kw...) function rebuild(trait::GI.AbstractTrait, geom, child_geoms; crs=GI.crs(geom), extent=nothing) T = GI.geointerface_geomtype(trait) - if GI.is3d(geom) - # The Boolean type parameters here indicate "3d-ness" and "measure" coordinate, respectively. - return T{true,false}(child_geoms; crs, extent) - else - return T{false,false}(child_geoms; crs, extent) - end + # Check the dimensionality of the first child geometry, since it may have changed + # NOTE that without this, 2D to 3D conversions will fail + hasZ = GI.is3d(first(child_geoms)) + hasM = GI.ismeasured(first(child_geoms)) + + return T{hasZ,hasM}(child_geoms; crs, extent) end diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index e0cd484bd..44897ff3b 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -70,6 +70,7 @@ include("transformations/segmentize.jl") include("transformations/simplify.jl") include("transformations/tuples.jl") include("transformations/transform.jl") +include("transformations/forcedims.jl") include("transformations/correction/geometry_correction.jl") include("transformations/correction/closed_ring.jl") include("transformations/correction/intersecting_polygons.jl") diff --git a/src/transformations/forcedims.jl b/src/transformations/forcedims.jl new file mode 100644 index 000000000..d614273b2 --- /dev/null +++ b/src/transformations/forcedims.jl @@ -0,0 +1,36 @@ +#= +# Force dimensions (xy, xyz) + +These functions force the geometry to be 2D or 3D. They work on any geometry, vector of geometries, feature collection, or table! + +They're implemented by `apply` pretty simply. +=# + +export forcexy, forcexyz + +""" + forcexy(geom) + +Force the geometry to be 2D. Works on any geometry, vector of geometries, feature collection, or table! +""" +function forcexy(geom) + return apply(GI.PointTrait(), geom) do point + (GI.x(point), GI.y(point)) + end +end + +""" + forcexyz(geom, z = 0) + +Force the geometry to be 3D. Works on any geometry, vector of geometries, feature collection, or table! + +The `z` parameter is the default z value - if a point has no z value, it will be set to this value. +If it does, then the z value will be kept. +""" +function forcexyz(geom, z = 0) + return apply(GI.PointTrait(), geom) do point + x, y = GI.x(point), GI.y(point) + z = GI.is3d(geom) ? GI.z(point) : z + (x, y, z) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index d4d25637a..3fc144484 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -27,6 +27,7 @@ include("helpers.jl") @safetestset "Simplify" begin include("transformations/simplify.jl") end @safetestset "Segmentize" begin include("transformations/segmentize.jl") end @safetestset "Transform" begin include("transformations/transform.jl") end +@safetestset "Force Dimensions" begin include("transformations/forcedims.jl") end @safetestset "Geometry Correction" begin include("transformations/correction/geometry_correction.jl") end @safetestset "Closed Rings" begin include("transformations/correction/closed_ring.jl") end @safetestset "Intersecting Polygons" begin include("transformations/correction/intersecting_polygons.jl") end diff --git a/test/transformations/forcedims.jl b/test/transformations/forcedims.jl new file mode 100644 index 000000000..61e4f9e3e --- /dev/null +++ b/test/transformations/forcedims.jl @@ -0,0 +1,35 @@ +import GeometryOps as GO +import GeoInterface as GI +using ..TestHelpers + +using Test + +geom = GI.Polygon([GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)]), +GI.LinearRing([(3, 4), (5, 6), (6, 7), (3, 4)])]) + +@testset_implementations "force dimensions" begin + + # Test forcexy on 3D geometry + geom3d = GO.transform($geom) do p + (GI.x(p), GI.y(p), 1.0) + end + @test GI.is3d(geom3d) + geom2d = GO.forcexy(geom3d) + @test !GI.is3d(geom2d) + @test GO.equals(geom2d, geom) + + # Test forcexyz with default z + geom3d_default = GO.forcexyz($geom) + @test GI.is3d(geom3d_default) + points3d = collect(GO.flatten(GI.PointTrait, geom3d_default)) + @test all(p -> GI.z(p) == 0, points3d) + + # Test forcexyz with custom z + geom3d_custom = GO.forcexyz($geom, 5.0) + @test GI.is3d(geom3d_custom) + points3d_custom = collect(GO.flatten(GI.PointTrait, geom3d_custom)) + @test all(p -> GI.z(p) == 5.0, points3d_custom) + + # Test forcexyz preserves existing z values + @test GO.equals(GO.forcexyz(geom3d), geom3d) +end diff --git a/test/transformations/transform.jl b/test/transformations/transform.jl index 40a700e65..770e39786 100644 --- a/test/transformations/transform.jl +++ b/test/transformations/transform.jl @@ -13,3 +13,17 @@ geom = GI.Polygon([GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)]), f = CoordinateTransformations.Translation(3.5, 1.5) @test GO.transform(f, $geom) == translated end + +@testset_implementations "transform 2D to 3D" begin + flat_points_raw = collect(GO.flatten(GI.PointTrait, $geom)) + flat_points_transformed = map(flat_points_raw) do p + (GI.x(p), GI.y(p), hypot(GI.x(p), GI.y(p))) + end + + geom_transformed = GO.transform($geom) do p + (GI.x(p), GI.y(p), hypot(GI.x(p), GI.y(p))) + end + @test collect(GO.flatten(GI.PointTrait, geom_transformed)) == flat_points_transformed + @test GI.is3d(geom_transformed) + @test !GI.ismeasured(geom_transformed) +end From c70d0071cb468dd45dedc1ecaed51efa9471a8d3 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Fri, 10 Jan 2025 14:40:28 -0100 Subject: [PATCH 26/30] Bugfix the type of `z` in forcexyz --- src/transformations/forcedims.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/transformations/forcedims.jl b/src/transformations/forcedims.jl index d614273b2..10340d1cc 100644 --- a/src/transformations/forcedims.jl +++ b/src/transformations/forcedims.jl @@ -30,7 +30,7 @@ If it does, then the z value will be kept. function forcexyz(geom, z = 0) return apply(GI.PointTrait(), geom) do point x, y = GI.x(point), GI.y(point) - z = GI.is3d(geom) ? GI.z(point) : z + z = GI.is3d(geom) ? GI.z(point) : convert(typeof(x), z) (x, y, z) end end From 1663f7082a56c272dc6453048823dbe032232674 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Fri, 10 Jan 2025 14:40:42 -0100 Subject: [PATCH 27/30] Declare compatibility with GeometryBasics v0.5 as well --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 1be7a8cb8..a4d4dc4d6 100644 --- a/Project.toml +++ b/Project.toml @@ -33,7 +33,7 @@ DelaunayTriangulation = "1.0.4" ExactPredicates = "2.2.8" FlexiJoins = "0.1.30" GeoInterface = "1.2" -GeometryBasics = "0.4.7" +GeometryBasics = "0.4.7, 0.5" GeometryOpsCore = "=0.1.2" LibGEOS = "0.9.2" LinearAlgebra = "1" From f409149a641433b167a0b251191cf439183a51fe Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Sun, 12 Jan 2025 07:42:04 -0500 Subject: [PATCH 28/30] Fix doctest failure when showing GeoInterface geometry --- src/utils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/utils.jl b/src/utils.jl index 82a1e664d..540a7211e 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -40,7 +40,7 @@ import GeometryOps as GO, GeoInterface as GI poly = GI.Polygon([[(-2.275543, 53.464547), (-2.275543, 53.489271), (-2.215118, 53.489271), (-2.215118, 53.464547), (-2.275543, 53.464547)]]) GO.polygon_to_line(poly) # output -GeoInterface.Wrappers.LineString{false, false, Vector{Tuple{Float64, Float64}}, Nothing, Nothing}([(-2.275543, 53.464547), (-2.275543, 53.489271), (-2.215118, 53.489271), (-2.215118, 53.464547), (-2.275543, 53.464547)], nothing, nothing) +GeoInterface.Wrappers.LineString{false, false}([(-2.275543, 53.464547), … (3) … , (-2.275543, 53.464547)]) ``` """ function polygon_to_line(poly) From 1f8ff9445af937cb4338d1c0b164b0254da3c8aa Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 13 Jan 2025 04:20:56 -0500 Subject: [PATCH 29/30] Add a brief writeup on what manifolds are to the docs --- docs/make.jl | 1 + docs/src/explanations/manifolds.md | 52 ++++++++++++++++++++++++++++++ 2 files changed, 53 insertions(+) create mode 100644 docs/src/explanations/manifolds.md diff --git a/docs/make.jl b/docs/make.jl index fb2ac0baf..a7ce8a1a5 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -113,6 +113,7 @@ makedocs(; "Explanations" => [ "Paradigms" => "explanations/paradigms.md", "Peculiarities" => "explanations/peculiarities.md", + "Manifolds" => "explanations/manifolds.md", ], "Source code" => literate_pages, ], diff --git a/docs/src/explanations/manifolds.md b/docs/src/explanations/manifolds.md new file mode 100644 index 000000000..fd0ffb336 --- /dev/null +++ b/docs/src/explanations/manifolds.md @@ -0,0 +1,52 @@ +# Manifolds + +A manifold is, mathematically, a description of some space that is locally Euclidean (i.e., locally flat). +All geographic projections, and the surface of the sphere and ellipsoid, fall under this category of space - +and these are all the spaces that are relevant to geographic geometry. + +## What manifolds are available? + +GeometryOps has three [`Manifold`](@ref) types: [`Planar`](@ref), [`Spherical`](@ref), and [`Geodesic`](@ref). + +- `Planar()` is, as the name suggests, a perfectly Cartesian, usually 2-dimensional, space. The shortest path from one point to another is a straight line. +- `Spherical(; radius)` describes points on the surface of a sphere of a given radius. + The most convenient sphere for geometry processing is the unit sphere, but one can also use + the sphere of the Earth for e.g. projections. +- `Geodesic(; semimajor_axis, inv_flattening)` describes points on the surface of a flattened ellipsoid, + similar to the Earth. The parameters describe the curvature and shape of the ellipsoid, and are equivalent + to the flags `+a` and `+f` in Proj's ellipsoid specification. The default values are the values of the WGS84 + ellipsoid. + + For `Geodesic`, we need an `AbstractGeodesic` that can wrap representations from Proj.jl and SphericalGeodesics.jl. + +The idea here is that the manifold describes how the geometry needs to be treated. + +## Why this is needed + +The classical problem this is intended to solve is that in GIS, latitude and longitude coordinates +are often treated as planar coordinates, when they in fact live on the sphere/ellipsoid, and must be +treated as such. For example, computing the area of the USA on the lat/long plane yields a result of `1116`, +which is plainly nonsensical. + +## How this is done + +In order to avoid this, we've introduced three complementary CRS-related systems to the JuliaGeo ecosystem. + +1. GeoInterface's `crstrait`. This is a method that returns the ideal CRS _type_ of a geometry, either Cartesian or Geographic. +2. Proj's `PreparedCRS` type, which extracts ellipsoid parameters and the nature of the projection from a coordinate reference system, and + caches the results in a struct. This allows GeometryOps to quickly determine the correct manifold to use for a given geometry. +3. GeometryOps's `Manifold` type, which defines the surface on which to perform operations. This is what allows GeometryOps to perform + calculations correctly depending on the nature of the geometry. + + +The way this flow works, is that when you load a geometry using GeoDataFrames, its CRS is extracted and parsed into a `PreparedCRS` type. +This is then used to determine the manifold to use for the geometry, and the geometry is converted to the manifold's coordinate system. + +There is a table of known geographic coordinate systems in GeoFormatTypes.jl, and anything else is assumed to be +a Cartesian or planar coordinate system. CRStrait is used as the cheap determinant, but PreparedCRS is more general and better to use if possible. + +When GeometryOps sees a geometry, it first checks its CRS to see if it is a geographic coordinate system. If it is, it uses the `PreparedCRS`, or falls back to `crstrait` and geographic defaults to determine the manifold. + +## Algorithms and manifolds + +Algorithms define what operation is performed on the geometry; however, the choice of algorithm can also depend on the manifold. L'Huilier's algorithm for the area of a polygon is not applicable to the plane, but is applicable to either the sphere or ellipsoid, for example. \ No newline at end of file From d0c0806e79c0bb293b5e6bacc2520b12f93a6ed1 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 13 Jan 2025 04:26:16 -0500 Subject: [PATCH 30/30] Use the GeoInterface show-fix branch --- .github/workflows/CI.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index cbbb3e165..ff1e526d6 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -60,7 +60,7 @@ jobs: with: version: '1' - name: Build and add versions - run: julia --project=docs -e 'using Pkg; Pkg.develop([PackageSpec(path = "."), PackageSpec(path = joinpath(".", "GeometryOpsCore"))]); Pkg.add([PackageSpec(name = "GeoMakie", rev = "master")])' + run: julia --project=docs -e 'using Pkg; Pkg.develop([PackageSpec(path = "."), PackageSpec(path = joinpath(".", "GeometryOpsCore"))]); Pkg.add([PackageSpec(name = "GeoMakie", rev = "master"), PackageSpec(name="GeoInterface", rev="bugfix_vars")])' - uses: julia-actions/julia-docdeploy@v1 with: install-package: false