[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"), PackageSpec(name="GeoInterface", rev="bugfix_vars")])'
       - 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/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

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

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"
@@ -32,7 +33,8 @@ 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"
 Proj = "1"

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,
     ],

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.

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

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
+function symdifference end

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) : convert(typeof(x), z)
+        (x, y, z)
+    end
+end

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

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)