diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 84a14a375..cbbb3e165 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: @@ -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.develop(; path = 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.develop([PackageSpec(path = "."), PackageSpec(path = joinpath(".", "GeometryOpsCore"))]); Pkg.add([PackageSpec(name = "GeoMakie", rev = "master")])' - uses: julia-actions/julia-docdeploy@v1 with: install-package: false 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 diff --git a/Project.toml b/Project.toml index 4bb1315b7..1be7a8cb8 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" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" SortTileRecursiveTree = "746ee33f-1797-42c2-866d-db2fce69d14d" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" @@ -33,6 +34,7 @@ ExactPredicates = "2.2.8" FlexiJoins = "0.1.30" GeoInterface = "1.2" GeometryBasics = "0.4.7" +GeometryOpsCore = "=0.1.2" LibGEOS = "0.9.2" LinearAlgebra = "1" Proj = "1" 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 diff --git a/ext/GeometryOpsProjExt/segmentize.jl b/ext/GeometryOpsProjExt/segmentize.jl index 308c97188..a7783ea02 100644 --- a/ext/GeometryOpsProjExt/segmentize.jl +++ b/ext/GeometryOpsProjExt/segmentize.jl @@ -1,20 +1,37 @@ # 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)) 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, ::Union{GI.LineStringTrait, GI.LinearRingTrait}; max_distance) + 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}[] + 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/GeometryOps.jl b/src/GeometryOps.jl index f524fd87d..e0cd484bd 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -2,16 +2,16 @@ module GeometryOps -include("../GeometryOpsCore/src/GeometryOpsCore.jl") # TODO: replace this with `using GeometryOpsCore` -import .GeometryOpsCore -for name in setdiff(names(GeometryOpsCore, all = true), (:eval, :var"#eval", :include, :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 +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 TraitTarget, Manifold, Planar, Spherical, Geodesic, apply, applyreduce, flatten, reconstruct, rebuild, unwrap using GeoInterface using GeometryBasics 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 diff --git a/src/transformations/segmentize.jl b/src/transformations/segmentize.jl index 623637bc5..df5882406 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,25 +156,38 @@ 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(Planar(), geom; max_distance, threaded = _booltype(threaded)) end + +# allow three-arg method as well, just in case +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()) + @assert max_distance > 0 "`max_distance` should be positive and nonzero! Found $(method.max_distance)." + _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()) + @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)." - segmentize_function = Base.Fix1(_segmentize, method) - return apply(segmentize_function, TraitTarget(GI.LinearRingTrait(), GI.LineStringTrait()), geom; threaded) + new_method = method isa LinearSegments ? Planar() : Geodesic() + segmentize(new_method, geom; max_distance = method.max_distance, threaded) end _segmentize(method, geom) = _segmentize(method, geom, GI.trait(geom)) @@ -182,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{LinearSegments, GeodesicSegments}, geom, T::Union{GI.LineStringTrait, GI.LinearRingTrait}) +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}[] @@ -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!(::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 > 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..10d07c7f3 100644 --- a/test/transformations/segmentize.jl +++ b/test/transformations/segmentize.jl @@ -31,39 +31,38 @@ 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 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.LinearSegments(; max_distance), $lr) + 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.GeodesicSegments(; max_distance), $lr) + @test_nowarn segmentized = GO.segmentize(GO.Geodesic(), $lr; max_distance) end end