[WIP] Use GeometryOpsCore for real

.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

GeometryOpsCore/Project.toml

@@ -1,7 +1,7 @@
 name = "GeometryOpsCore"
 uuid = "05efe853-fabf-41c8-927e-7063c8b9f013"
 authors = ["Anshul Singhvi <[email protected]>", "Rafael Schouten <[email protected]>", "Skylar Gering <[email protected]>", "and contributors"]
-version = "0.1.1"
+version = "0.1.2"
 
 [deps]
 DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"

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
 

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"

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

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))

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

src/transformations/segmentize.jl

@@ -146,7 +146,7 @@
 # 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,51 +156,64 @@
 # ## 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))
 #= 
 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}[]
     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)
+        _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))

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