Use GeometryOpsCore for real (#223)

Author: asinghvi17

.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

@@ -126,8 +126,8 @@ end
 @inline function apply(
     f::F, target, geom; calc_extent=false, threaded=false, kw...
 ) where F
-    threaded = _booltype(threaded)
-    calc_extent = _booltype(calc_extent)
+    threaded = booltype(threaded)
+    calc_extent = booltype(calc_extent)
     _apply(f, TraitTarget(target), geom; threaded, calc_extent, kw...)
 end
 
@@ -139,7 +139,7 @@ end
     # For an Array there is nothing else to do but map `_apply` over all values
     # _maptasks may run this level threaded if `threaded==true`,
     # but deeper `_apply` called in the closure will not be threaded
-    apply_to_array(i) = _apply(f, target, A[i]; threaded=_False(), kw...)
+    apply_to_array(i) = _apply(f, target, A[i]; threaded=False(), kw...)
     _maptasks(apply_to_array, eachindex(A), threaded)
 end
 # There is no trait and this is not an AbstractArray.
@@ -150,7 +150,7 @@ end
     if Tables.istable(iterable)
     _apply_table(f, target, iterable; threaded, kw...)
     else # this is probably some form of iterable...
-        if threaded isa _True
+        if threaded isa True
             # `collect` first so we can use threads
             _apply(f, target, collect(iterable); threaded, kw...)
         else
@@ -226,14 +226,14 @@ end
 # Rewrap all FeatureCollectionTrait feature collections as GI.FeatureCollection
 # Maybe use threads to call _apply on component features
 @inline function _apply(f::F, target, ::GI.FeatureCollectionTrait, fc;
-    crs=GI.crs(fc), calc_extent=_False(), threaded
+    crs=GI.crs(fc), calc_extent=False(), threaded
 ) where F
 
     # Run _apply on all `features` in the feature collection, possibly threaded
     apply_to_feature(i) =
-        _apply(f, target, GI.getfeature(fc, i); crs, calc_extent, threaded=_False())::GI.Feature
+        _apply(f, target, GI.getfeature(fc, i); crs, calc_extent, threaded=False())::GI.Feature
     features = _maptasks(apply_to_feature, 1:GI.nfeature(fc), threaded)
-    if calc_extent isa _True
+    if calc_extent isa True
         # Calculate the extent of the features
         extent = mapreduce(GI.extent, Extents.union, features)
         # Return a FeatureCollection with features, crs and calculated extent
@@ -245,13 +245,13 @@ end
 end
 # Rewrap all FeatureTrait features as GI.Feature, keeping the properties
 @inline function _apply(f::F, target, ::GI.FeatureTrait, feature;
-    crs=GI.crs(feature), calc_extent=_False(), threaded
+    crs=GI.crs(feature), calc_extent=False(), threaded
 ) where F
     # Run _apply on the contained geometry
     geometry = _apply(f, target, GI.geometry(feature); crs, calc_extent, threaded)
     # Get the feature properties
     properties = GI.properties(feature)
-    if calc_extent isa _True
+    if calc_extent isa True
         # Calculate the extent of the geometry
         extent = GI.extent(geometry)
         # Return a new Feature with the new geometry and calculated extent, but the original properties and crs
@@ -264,36 +264,36 @@ end
 # Reconstruct nested geometries,
 # maybe using threads to call _apply on component geoms
 @inline function _apply(f::F, target, trait, geom;
-    crs=GI.crs(geom), calc_extent=_False(), threaded
+    crs=GI.crs(geom), calc_extent=False(), threaded
 )::(GI.geointerface_geomtype(trait)) where F
     # Map `_apply` over all sub geometries of `geom`
     # to create a new vector of geometries
     # TODO handle zero length
-    apply_to_geom(i) = _apply(f, target, GI.getgeom(geom, i); crs, calc_extent, threaded=_False())
+    apply_to_geom(i) = _apply(f, target, GI.getgeom(geom, i); crs, calc_extent, threaded=False())
     geoms = _maptasks(apply_to_geom, 1:GI.ngeom(geom), threaded)
     return _apply_inner(geom, geoms, crs, calc_extent)
 end
 @inline function _apply(f::F, target::TraitTarget{<:PointTrait}, trait::GI.PolygonTrait, geom;
-    crs=GI.crs(geom), calc_extent=_False(), threaded
+    crs=GI.crs(geom), calc_extent=False(), threaded
 )::(GI.geointerface_geomtype(trait)) where F
     # We need to force rebuilding a LinearRing not a LineString
     geoms = _maptasks(1:GI.ngeom(geom), threaded) do i
         lr = GI.getgeom(geom, i)
         points = map(GI.getgeom(lr)) do p
-            _apply(f, target, p; crs, calc_extent, threaded=_False())
+            _apply(f, target, p; crs, calc_extent, threaded=False())
         end
         _linearring(_apply_inner(lr, points, crs, calc_extent))
     end
     return _apply_inner(geom, geoms, crs, calc_extent)
 end
-function _apply_inner(geom, geoms, crs, calc_extent::_True)
+function _apply_inner(geom, geoms, crs, calc_extent::True)
     # Calculate the extent of the sub geometries
     extent = mapreduce(GI.extent, Extents.union, geoms)
     # Return a new geometry of the same trait as `geom`,
     # holding the new `geoms` with `crs` and calculated extent
     return rebuild(geom, geoms; crs, extent)
 end
-function _apply_inner(geom, geoms, crs, calc_extent::_False)
+function _apply_inner(geom, geoms, crs, calc_extent::False)
     # Return a new geometry of the same trait as `geom`, holding the new `geoms` with `crs`
     return rebuild(geom, geoms; crs)
 end
@@ -322,7 +322,7 @@ using Base.Threads: nthreads, @threads, @spawn
 # Threading utility, modified Mason Protters threading PSA
 # run `f` over ntasks, where f receives an AbstractArray/range
 # of linear indices
-@inline function _maptasks(f::F, taskrange, threaded::_True)::Vector where F
+@inline function _maptasks(f::F, taskrange, threaded::True)::Vector where F
     ntasks = length(taskrange)
     # Customize this as needed.
     # More tasks have more overhead, but better load balancing
@@ -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/applyreduce.jl

@@ -34,24 +34,24 @@ feature collections and nested geometries.
 @inline function applyreduce(
     f::F, op::O, target, geom; threaded=false, init=nothing
 ) where {F, O}
-    threaded = _booltype(threaded)
+    threaded = booltype(threaded)
     _applyreduce(f, op, TraitTarget(target), geom; threaded, init)
 end
 
 @inline _applyreduce(f::F, op::O, target, geom; threaded, init) where {F, O} =
     _applyreduce(f, op, target, GI.trait(geom), geom; threaded, init)
 # Maybe use threads reducing over arrays
 @inline function _applyreduce(f::F, op::O, target, ::Nothing, A::AbstractArray; threaded, init) where {F, O}
-    applyreduce_array(i) = _applyreduce(f, op, target, A[i]; threaded=_False(), init)
+    applyreduce_array(i) = _applyreduce(f, op, target, A[i]; threaded=False(), init)
     _mapreducetasks(applyreduce_array, op, eachindex(A), threaded; init)
 end
 # Try to applyreduce over iterables
 @inline function _applyreduce(f::F, op::O, target, ::Nothing, iterable::IterableType; threaded, init) where {F, O, IterableType}
     if Tables.istable(iterable)
         _applyreduce_table(f, op, target, iterable; threaded, init)
     else
-        applyreduce_iterable(i) = _applyreduce(f, op, target, i; threaded=_False(), init)
-        if threaded isa _True # Try to `collect` and reduce over the vector with threads
+        applyreduce_iterable(i) = _applyreduce(f, op, target, i; threaded=False(), init)
+        if threaded isa True # Try to `collect` and reduce over the vector with threads
             _applyreduce(f, op, target, collect(iterable); threaded, init)
         else
             # Try to `mapreduce` the iterable as-is
@@ -77,23 +77,23 @@ function _applyreduce_table(f::F, op::O, target::GI.FeatureTrait, iterable::Iter
 end
 # Maybe use threads reducing over features of feature collections
 @inline function _applyreduce(f::F, op::O, target, ::GI.FeatureCollectionTrait, fc; threaded, init) where {F, O}
-    applyreduce_fc(i) = _applyreduce(f, op, target, GI.getfeature(fc, i); threaded=_False(), init)
+    applyreduce_fc(i) = _applyreduce(f, op, target, GI.getfeature(fc, i); threaded=False(), init)
     _mapreducetasks(applyreduce_fc, op, 1:GI.nfeature(fc), threaded; init)
 end
 # Features just applyreduce to their geometry
 @inline _applyreduce(f::F, op::O, target, ::GI.FeatureTrait, feature; threaded, init) where {F, O} =
     _applyreduce(f, op, target, GI.geometry(feature); threaded, init)
 # Maybe use threads over components of nested geometries
 @inline function _applyreduce(f::F, op::O, target, trait, geom; threaded, init) where {F, O}
-    applyreduce_geom(i) = _applyreduce(f, op, target, GI.getgeom(geom, i); threaded=_False(), init)
+    applyreduce_geom(i) = _applyreduce(f, op, target, GI.getgeom(geom, i); threaded=False(), init)
     _mapreducetasks(applyreduce_geom, op, 1:GI.ngeom(geom), threaded; init)
 end
 # Don't thread over points it won't pay off
 @inline function _applyreduce(
     f::F, op::O, target, trait::Union{GI.LinearRing,GI.LineString,GI.MultiPoint}, geom;
     threaded, init
 ) where {F, O}
-    _applyreduce(f, op, target, GI.getgeom(geom); threaded=_False(), init)
+    _applyreduce(f, op, target, GI.getgeom(geom); threaded=False(), init)
 end
 # Apply f to the target
 @inline function _applyreduce(f::F, op::O, ::TraitTarget{Target}, ::Trait, x; kw...) where {F,O,Target,Trait<:Target} 
@@ -124,7 +124,7 @@ import Base.Threads: nthreads, @threads, @spawn
 # 
 # If you absolutely need a single chunk, then `threaded = false` will always decompose
 # to straight `mapreduce` without grouping.
-@inline function _mapreducetasks(f::F, op, taskrange, threaded::_True; init) where F
+@inline function _mapreducetasks(f::F, op, taskrange, threaded::True; init) where F
     ntasks = length(taskrange)
     # Customize this as needed.
     # More tasks have more overhead, but better load balancing
@@ -144,6 +144,6 @@ import Base.Threads: nthreads, @threads, @spawn
     # Finally we join the results into a new vector
     return mapreduce(fetch, op, tasks; init)
 end
-Base.@assume_effects :foldable function _mapreducetasks(f::F, op, taskrange, threaded::_False; init) where F
+Base.@assume_effects :foldable function _mapreducetasks(f::F, op, taskrange, threaded::False; init) where F
     mapreduce(f, op, taskrange; init)
 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

@@ -22,7 +22,7 @@ Currently, those circumstances are `area` and `segmentize`, but this could be ex
 
 export Planar, Spherical, Geodesic
 export TraitTarget
-export BoolsAsTypes, _True, _False, _booltype
+export BoolsAsTypes, True, False, booltype
 
 """
     abstract type Manifold
@@ -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
 
@@ -128,7 +130,7 @@ the compiler to separate threaded and non-threaded code paths.
 Note that if we didn't include the parent abstract type, this would have been really 
 type unstable, since the compiler couldn't tell what would be returned!
 
-We had to add the type annotation on the `_booltype(::Bool)` method for this reason as well.
+We had to add the type annotation on the `booltype(::Bool)` method for this reason as well.
 
 TODO: should we switch to `Static.jl`?
 =#
@@ -140,25 +142,25 @@ TODO: should we switch to `Static.jl`?
 abstract type BoolsAsTypes end
 
 """
-    struct _True <: BoolsAsTypes
+    struct True <: BoolsAsTypes
 
 A struct that means `true`.
 """
-struct _True <: BoolsAsTypes end
+struct True <: BoolsAsTypes end
 
 """
-    struct _False <: BoolsAsTypes
+    struct False <: BoolsAsTypes
 
 A struct that means `false`.
 """
-struct _False <: BoolsAsTypes end
+struct False <: BoolsAsTypes end
 
 """
-    _booltype(x)
+    booltype(x)
 
 Returns a [`BoolsAsTypes`](@ref) from `x`, whether it's a boolean or a BoolsAsTypes.
 """
-function _booltype end
+function booltype end
 
-@inline _booltype(x::Bool)::BoolsAsTypes = x ? _True() : _False()
-@inline _booltype(x::BoolsAsTypes)::BoolsAsTypes = x
+@inline booltype(x::Bool)::BoolsAsTypes = x ? True() : False()
+@inline booltype(x::BoolsAsTypes)::BoolsAsTypes = x

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