From 0a2d50e4715974c74c76995ecd66a508a6fe3f67 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 2 Oct 2024 15:18:43 -0700 Subject: [PATCH 01/10] Factor all core stuff out into a GeometryOpsCore package This basically guts the `types.jl` and `primitives.jl` files from GeometryOps and puts them into a new, low-dependency module that other packages can utilize for a convenient structure. --- .gitignore | 1 + GeometryOpsCore/Project.toml | 4 + GeometryOpsCore/README.md | 0 GeometryOpsCore/src/GeometryOpsCore.jl | 17 + GeometryOpsCore/src/apply.jl | 351 +++++++++++++ GeometryOpsCore/src/applyreduce.jl | 149 ++++++ GeometryOpsCore/src/keyword_docs.jl | 15 + GeometryOpsCore/src/other_primitives.jl | 172 +++++++ GeometryOpsCore/src/types.jl | 95 ++++ src/GeometryOps.jl | 11 + src/primitives.jl | 652 ------------------------ src/types.jl | 61 +-- 12 files changed, 816 insertions(+), 712 deletions(-) create mode 100644 GeometryOpsCore/Project.toml create mode 100644 GeometryOpsCore/README.md create mode 100644 GeometryOpsCore/src/GeometryOpsCore.jl create mode 100644 GeometryOpsCore/src/apply.jl create mode 100644 GeometryOpsCore/src/applyreduce.jl create mode 100644 GeometryOpsCore/src/keyword_docs.jl create mode 100644 GeometryOpsCore/src/other_primitives.jl create mode 100644 GeometryOpsCore/src/types.jl diff --git a/.gitignore b/.gitignore index f0a2dc8bb..37d130b1a 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ /Manifest.toml +/GeometryOpsCore/Manifest.toml /docs/build/ /docs/src/source/ .vscode/ diff --git a/GeometryOpsCore/Project.toml b/GeometryOpsCore/Project.toml new file mode 100644 index 000000000..2b9dcca14 --- /dev/null +++ b/GeometryOpsCore/Project.toml @@ -0,0 +1,4 @@ +[deps] +DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" +GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f" +Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c" diff --git a/GeometryOpsCore/README.md b/GeometryOpsCore/README.md new file mode 100644 index 000000000..e69de29bb diff --git a/GeometryOpsCore/src/GeometryOpsCore.jl b/GeometryOpsCore/src/GeometryOpsCore.jl new file mode 100644 index 000000000..4084f13db --- /dev/null +++ b/GeometryOpsCore/src/GeometryOpsCore.jl @@ -0,0 +1,17 @@ +module GeometryOpsCore + +using Base.Threads: nthreads, @threads, @spawn + +import GeoInterface as GI + +using Tables +using DataAPI + +include("keyword_docs.jl") +include("types.jl") + +include("apply.jl") +include("applyreduce.jl") +include("other_primitives.jl") + +end \ No newline at end of file diff --git a/GeometryOpsCore/src/apply.jl b/GeometryOpsCore/src/apply.jl new file mode 100644 index 000000000..d7713195d --- /dev/null +++ b/GeometryOpsCore/src/apply.jl @@ -0,0 +1,351 @@ +# # `apply` + +export apply + +#= + +This file mainly defines the [`apply`](@ref) function. + +In general, the idea behind the `apply` framework is to take +as input any geometry, vector of geometries, or feature collection, +deconstruct it to the given trait target (any arbitrary GI.AbstractTrait +or `TraitTarget` union thereof, like `PointTrait` or `PolygonTrait`) +and perform some operation on it. Then, the geometry or structure is rebuilt. + +This allows for a simple and consistent framework within which users can +define their own operations trivially easily, and removes a lot of the +complexity involved with handling complex geometry structures. + +For example, a simple way to flip the x and y coordinates of a geometry is: + +```julia +flipped_geom = GO.apply(GI.PointTrait(), geom) do p + (GI.y(p), GI.x(p)) +end +``` + +As simple as that. There's no need to implement your own decomposition because it's done for you. + +Functions like [`flip`](@ref), [`reproject`](@ref), [`transform`](@ref), even [`segmentize`](@ref) and [`simplify`](@ref) have been implemented +using the `apply` framework. Similarly, [`centroid`](@ref), [`area`](@ref) and [`distance`](@ref) have been implemented using the +[`applyreduce`](@ref) framework. + +## Docstrings + +### Functions + +```@docs; collapse=true, canonical=false +apply +applyreduce +``` + +=# + +#= +## What is `apply`? + +`apply` applies some function to every geometry matching the `Target` +GeoInterface trait, in some arbitrarily nested object made up of: +- `AbstractArray`s (we also try to iterate other non-GeoInteface compatible object) +- `FeatureCollectionTrait` objects +- `FeatureTrait` objects +- `AbstractGeometryTrait` objects + +`apply` recursively calls itself through these nested +layers until it reaches objects with the `Target` GeoInterface trait. When found `apply` applies the function `f`, and stops. + +The outer recursive functions then progressively rebuild the object +using GeoInterface objects matching the original traits. + +If `PointTrait` is found but it is not the `Target`, an error is thrown. +This likely means the object contains a different geometry trait to +the target, such as `MultiPointTrait` when `LineStringTrait` was specified. + +To handle this possibility it may be necessary to make `Target` a +`Union` of traits found at the same level of nesting, and define methods +of `f` to handle all cases. + +Be careful making a union across "levels" of nesting, e.g. +`Union{FeatureTrait,PolygonTrait}`, as `_apply` will just never reach +`PolygonTrait` when all the polygons are wrapped in a `FeatureTrait` object. + +## Embedding: + +`extent` and `crs` can be embedded in all geometries, features, and +feature collections as part of `apply`. Geometries deeper than `Target` +will of course not have new `extent` or `crs` embedded. + +- `calc_extent` signals to recalculate an `Extent` and embed it. +- `crs` will be embedded as-is + +## Threading + +Threading is used at the outermost level possible - over +an array, feature collection, or e.g. a MultiPolygonTrait where +each `PolygonTrait` sub-geometry may be calculated on a different thread. + +Currently, threading defaults to `false` for all objects, but can be turned on +by passing the keyword argument `threaded=true` to `apply`. +=# + +""" + apply(f, target::Union{TraitTarget, GI.AbstractTrait}, obj; kw...) + +Reconstruct a geometry, feature, feature collection, or nested vectors of +either using the function `f` on the `target` trait. + +`f(target_geom) => x` where `x` also has the `target` trait, or a trait that can +be substituted. For example, swapping `PolgonTrait` to `MultiPointTrait` will fail +if the outer object has `MultiPolygonTrait`, but should work if it has `FeatureTrait`. + +Objects "shallower" than the target trait are always completely rebuilt, like +a `Vector` of `FeatureCollectionTrait` of `FeatureTrait` when the target +has `PolygonTrait` and is held in the features. These will always be GeoInterface +geometries/feature/feature collections. But "deeper" objects may remain +unchanged or be whatever GeoInterface compatible objects `f` returns. + +The result is a functionally similar geometry with values depending on `f`. + +$APPLY_KEYWORDS + +## Example + +Flipped point the order in any feature or geometry, or iterables of either: + +```julia +import GeoInterface as GI +import GeometryOps as GO +geom = GI.Polygon([GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)]), + GI.LinearRing([(3, 4), (5, 6), (6, 7), (3, 4)])]) + +flipped_geom = GO.apply(GI.PointTrait, geom) do p + (GI.y(p), GI.x(p)) +end +``` +""" +@inline function apply( + f::F, target, geom; calc_extent=false, threaded=false, kw... +) where F + threaded = _booltype(threaded) + calc_extent = _booltype(calc_extent) + _apply(f, TraitTarget(target), geom; threaded, calc_extent, kw...) +end + +# Call _apply again with the trait of `geom` +@inline _apply(f::F, target, geom; kw...) where F = + _apply(f, target, GI.trait(geom), geom; kw...) +# There is no trait and this is an AbstractArray - so just iterate over it calling _apply on the contents +@inline function _apply(f::F, target, ::Nothing, A::AbstractArray; threaded, kw...) where F + # 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...) + _maptasks(apply_to_array, eachindex(A), threaded) +end +# There is no trait and this is not an AbstractArray. +# Try to call _apply over it. We can't use threading +# as we don't know if we can can index into it. So just `map`. +@inline function _apply(f::F, target, ::Nothing, iterable::IterableType; threaded, kw...) where {F, IterableType} + # Try the Tables.jl interface first + if Tables.istable(iterable) + _apply_table(f, target, iterable; threaded, kw...) + else # this is probably some form of iterable... + if threaded isa _True + # `collect` first so we can use threads + _apply(f, target, collect(iterable); threaded, kw...) + else + apply_to_iterable(x) = _apply(f, target, x; kw...) + map(apply_to_iterable, iterable) + end + end +end +#= +Doing this inline in `_apply` is _heavily_ type unstable, so it's best to separate this +by a function barrier. + +This function operates `apply` on the `geometry` column of the table, and returns a new table +with the same schema, but with the new geometry column. + +This new table may be of the same type as the old one iff `Tables.materializer` is defined for +that table. If not, then a `NamedTuple` is returned. +=# +function _apply_table(f::F, target, iterable::IterableType; threaded, kw...) where {F, IterableType} + _get_col_pair(colname) = colname => Tables.getcolumn(iterable, colname) + # We extract the geometry column and run `apply` on it. + geometry_column = first(GI.geometrycolumns(iterable)) + new_geometry = _apply(f, target, Tables.getcolumn(iterable, geometry_column); threaded, kw...) + # Then, we obtain the schema of the table, + old_schema = Tables.schema(iterable) + # filter the geometry column out, + new_names = filter(Base.Fix1(!==, geometry_column), old_schema.names) + # and try to rebuild the same table as the best type - either the original type of `iterable`, + # or a named tuple which is the default fallback. + result = Tables.materializer(iterable)( + merge( + NamedTuple{(geometry_column,), Base.Tuple{typeof(new_geometry)}}((new_geometry,)), + NamedTuple(Iterators.map(_get_col_pair, new_names)) + ) + ) + # Finally, we ensure that metadata is propagated correctly. + # This can only happen if the original table supports metadata reads, + # and the result supports metadata writes. + if DataAPI.metadatasupport(typeof(result)).write + # Copy over all metadata from the original table to the new table, + # if the original table supports metadata reading. + if DataAPI.metadatasupport(IterableType).read + for (key, (value, style)) in DataAPI.metadata(iterable; style = true) + # Default styles are not preserved on data transformation, so we must skip them! + style == :default && continue + # We assume that any other style is preserved. + DataAPI.metadata!(result, key, value; style) + end + end + # We don't usually care about the original table's metadata for GEOINTERFACE namespaced + # keys, so we should set the crs and geometrycolumns metadata if they are present. + # Ensure that `GEOINTERFACE:geometrycolumns` and `GEOINTERFACE:crs` are set! + mdk = DataAPI.metadatakeys(result) + # If the user has asked for geometry columns to persist, they would be here, + # so we don't need to set them. + if !("GEOINTERFACE:geometrycolumns" in mdk) + # If the geometry columns are not already set, we need to set them. + DataAPI.metadata!(result, "GEOINTERFACE:geometrycolumns", (geometry_column,); style = :default) + end + # Force reset CRS always, since you can pass `crs` to `apply`. + new_crs = if haskey(kw, :crs) + kw[:crs] + else + GI.crs(iterable) # this will automatically check `GEOINTERFACE:crs` unless the type has a specialized implementation. + end + + DataAPI.metadata!(result, "GEOINTERFACE:crs", new_crs; style = :default) + end + + return result +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 +) 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 + features = _maptasks(apply_to_feature, 1:GI.nfeature(fc), threaded) + 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 + return GI.FeatureCollection(features; crs, extent) + else + # Return a FeatureCollection with features and crs + return GI.FeatureCollection(features; crs) + 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 +) 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 + # 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 + return GI.Feature(geometry; properties, crs, extent) + else + # Return a new Feature with the new geometry, but the original properties and crs + return GI.Feature(geometry; properties, crs) + end +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 +)::(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()) + 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 +)::(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()) + 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) + # 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) + # Return a new geometry of the same trait as `geom`, holding the new `geoms` with `crs` + return rebuild(geom, geoms; crs) +end +# Fail loudly if we hit PointTrait without running `f` +# (after PointTrait there is no further to dig with `_apply`) +# @inline _apply(f, ::TraitTarget{Target}, trait::GI.PointTrait, geom; crs=nothing, kw...) where Target = + # throw(ArgumentError("target $Target not found, but reached a `PointTrait` leaf")) +# Finally, these short methods are the main purpose of `apply`. +# The `Trait` is a subtype of the `Target` (or identical to it) +# So the `Target` is found. We apply `f` to geom and return it to previous +# _apply calls to be wrapped with the outer geometries/feature/featurecollection/array. +_apply(f::F, ::TraitTarget{Target}, ::Trait, geom; crs=GI.crs(geom), kw...) where {F,Target,Trait<:Target} = f(geom) +# Define some specific cases of this match to avoid method ambiguity +for T in ( + GI.PointTrait, GI.LinearRing, GI.LineString, + GI.MultiPoint, GI.FeatureTrait, GI.FeatureCollectionTrait +) + @eval _apply(f::F, target::TraitTarget{<:$T}, trait::$T, x; kw...) where F = f(x) +end + + +### `_maptasks` - flexible, threaded `map` + +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 + ntasks = length(taskrange) + # Customize this as needed. + # More tasks have more overhead, but better load balancing + tasks_per_thread = 2 + chunk_size = max(1, ntasks ÷ (tasks_per_thread * nthreads())) + # partition the range into chunks + task_chunks = Iterators.partition(taskrange, chunk_size) + # Map over the chunks + tasks = map(task_chunks) do chunk + # Spawn a task to process this chunk + @spawn begin + # Where we map `f` over the chunk indices + map(f, chunk) + end + end + + # Finally we join the results into a new vector + return mapreduce(fetch, vcat, tasks) +end +#= +Here we 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! +=# +Base.@assume_effects :foldable @inline function _maptasks(f::F, taskrange, threaded::_False)::Vector where F + map(f, taskrange) +end \ No newline at end of file diff --git a/GeometryOpsCore/src/applyreduce.jl b/GeometryOpsCore/src/applyreduce.jl new file mode 100644 index 000000000..5e21c3a2d --- /dev/null +++ b/GeometryOpsCore/src/applyreduce.jl @@ -0,0 +1,149 @@ +#= +# `applyreduce` +=# + +export applyreduce + +#= +This file mainly defines the [`applyreduce`](@ref) function. + +This performs `apply`, but then reduces the result after flattening instead of rebuilding the geometry. + + +In general, the idea behind the `apply` framework is to take +as input any geometry, vector of geometries, or feature collection, +deconstruct it to the given trait target (any arbitrary GI.AbstractTrait +or `TraitTarget` union thereof, like `PointTrait` or `PolygonTrait`) +and perform some operation on it. + +[`centroid`](@ref), [`area`](@ref) and [`distance`](@ref) have been implemented using the +[`applyreduce`](@ref) framework. +=# + +""" + applyreduce(f, op, target::Union{TraitTarget, GI.AbstractTrait}, obj; threaded) + +Apply function `f` to all objects with the `target` trait, +and reduce the result with an `op` like `+`. + +The order and grouping of application of `op` is not guaranteed. + +If `threaded==true` threads will be used over arrays and iterables, +feature collections and nested geometries. +""" +@inline function applyreduce( + f::F, op::O, target, geom; threaded=false, init=nothing +) where {F, O} + 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) + _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(f, op, target, collect(iterable); threaded, init) + else + # Try to `mapreduce` the iterable as-is + mapreduce(applyreduce_iterable, op, iterable; init) + end + end +end +# In this case, we don't reconstruct the table, but only operate on the geometry column. +function _applyreduce_table(f::F, op::O, target, iterable::IterableType; threaded, init) where {F, O, IterableType} + # We extract the geometry column and run `applyreduce` on it. + geometry_column = first(GI.geometrycolumns(iterable)) + return _applyreduce(f, op, target, Tables.getcolumn(iterable, geometry_column); threaded, init) +end +# If `applyreduce` wants features, then applyreduce over the rows as `GI.Feature`s. +function _applyreduce_table(f::F, op::O, target::GI.FeatureTrait, iterable::IterableType; threaded, init) where {F, O, IterableType} + # We extract the geometry column and run `apply` on it. + geometry_column = first(GI.geometrycolumns(iterable)) + property_names = Iterators.filter(!=(geometry_column), Tables.schema(iterable).names) + features = map(Tables.rows(iterable)) do row + GI.Feature(Tables.getcolumn(row, geometry_column), properties=NamedTuple(Iterators.map(Base.Fix1(_get_col_pair, row), property_names))) + end + return _applyreduce(f, op, target, features; threaded, init) +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) + _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) + _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) +end +# Apply f to the target +@inline function _applyreduce(f::F, op::O, ::TraitTarget{Target}, ::Trait, x; kw...) where {F,O,Target,Trait<:Target} + f(x) +end +# Fail if we hit PointTrait +# _applyreduce(f, op, target::TraitTarget{Target}, trait::PointTrait, geom; kw...) where Target = + # throw(ArgumentError("target $target not found")) +# Specific cases to avoid method ambiguity +for T in ( + GI.PointTrait, GI.LinearRing, GI.LineString, + GI.MultiPoint, GI.FeatureTrait, GI.FeatureCollectionTrait +) + @eval _applyreduce(f::F, op::O, ::TraitTarget{<:$T}, trait::$T, x; kw...) where {F, O} = f(x) +end + +### `_mapreducetasks` - flexible, threaded mapreduce + +import Base: nthreads, @threads, @spawn + +# Threading utility, modified Mason Protters threading PSA +# run `f` over ntasks, where f receives an AbstractArray/range +# of linear indices +# +# WARNING: this will not work for mean/median - only ops +# where grouping is possible. That's because the implementation operates +# in chunks, and not globally. +# +# 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 + ntasks = length(taskrange) + # Customize this as needed. + # More tasks have more overhead, but better load balancing + tasks_per_thread = 2 + chunk_size = max(1, ntasks ÷ (tasks_per_thread * nthreads())) + # partition the range into chunks + task_chunks = Iterators.partition(taskrange, chunk_size) + # Map over the chunks + tasks = map(task_chunks) do chunk + # Spawn a task to process this chunk + @spawn begin + # Where we map `f` over the chunk indices + mapreduce(f, op, chunk; init) + end + end + + # 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 + mapreduce(f, op, taskrange; init) +end diff --git a/GeometryOpsCore/src/keyword_docs.jl b/GeometryOpsCore/src/keyword_docs.jl new file mode 100644 index 000000000..a4067d79b --- /dev/null +++ b/GeometryOpsCore/src/keyword_docs.jl @@ -0,0 +1,15 @@ +#= +# Keyword docs + +This file defines common keyword documentation, that can be spliced into docstrings. +=# + +const THREADED_KEYWORD = "- `threaded`: `true` or `false`. Whether to use multithreading. Defaults to `false`." +const CRS_KEYWORD = "- `crs`: The CRS to attach to geometries. Defaults to `nothing`." +const CALC_EXTENT_KEYWORD = "- `calc_extent`: `true` or `false`. Whether to calculate the extent. Defaults to `false`." + +const APPLY_KEYWORDS = """ +$THREADED_KEYWORD +$CRS_KEYWORD +$CALC_EXTENT_KEYWORD +""" \ No newline at end of file diff --git a/GeometryOpsCore/src/other_primitives.jl b/GeometryOpsCore/src/other_primitives.jl new file mode 100644 index 000000000..b0fd04de7 --- /dev/null +++ b/GeometryOpsCore/src/other_primitives.jl @@ -0,0 +1,172 @@ +#= +# Other primitives (unwrap, flatten, etc) + +This file defines the following primitives: +```@docs; collapsed=true, canonical=false +unwrap +flatten +reconstruct +rebuild +``` +=# + +""" + unwrap(target::Type{<:AbstractTrait}, obj) + unwrap(f, target::Type{<:AbstractTrait}, obj) + +Unwrap the object to vectors, down to the target trait. + +If `f` is passed in it will be applied to the target geometries +as they are found. +""" +function unwrap end +unwrap(target::Type, geom) = unwrap(identity, target, geom) +# Add dispatch argument for trait +unwrap(f, target::Type, geom) = unwrap(f, target, GI.trait(geom), geom) +# Try to unwrap over iterables +unwrap(f, target::Type, ::Nothing, iterable) = + map(x -> unwrap(f, target, x), iterable) +# Rewrap feature collections +unwrap(f, target::Type, ::GI.FeatureCollectionTrait, fc) = + map(x -> unwrap(f, target, x), GI.getfeature(fc)) +unwrap(f, target::Type, ::GI.FeatureTrait, feature) = + unwrap(f, target, GI.geometry(feature)) +unwrap(f, target::Type, trait, geom) = map(g -> unwrap(f, target, g), GI.getgeom(geom)) +# Apply f to the target geometry +unwrap(f, ::Type{Target}, ::Trait, geom) where {Target,Trait<:Target} = f(geom) +# Fail if we hit PointTrait +unwrap(f, target::Type, trait::GI.PointTrait, geom) = + throw(ArgumentError("target $target not found, but reached a `PointTrait` leaf")) +# Specific cases to avoid method ambiguity +unwrap(f, target::Type{GI.PointTrait}, trait::GI.PointTrait, geom) = f(geom) +unwrap(f, target::Type{GI.FeatureTrait}, ::GI.FeatureTrait, feature) = f(feature) +unwrap(f, target::Type{GI.FeatureCollectionTrait}, ::GI.FeatureCollectionTrait, fc) = f(fc) + +""" + flatten(target::Type{<:GI.AbstractTrait}, obj) + flatten(f, target::Type{<:GI.AbstractTrait}, obj) + +Lazily flatten any `AbstractArray`, iterator, `FeatureCollectionTrait`, +`FeatureTrait` or `AbstractGeometryTrait` object `obj`, so that objects +with the `target` trait are returned by the iterator. + +If `f` is passed in it will be applied to the target geometries. +""" +flatten(::Type{Target}, geom) where {Target<:GI.AbstractTrait} = flatten(identity, Target, geom) +flatten(f, ::Type{Target}, geom) where {Target<:GI.AbstractTrait} = _flatten(f, Target, geom) + +_flatten(f, ::Type{Target}, geom) where Target = _flatten(f, Target, GI.trait(geom), geom) +# Try to flatten over iterables +function _flatten(f, ::Type{Target}, ::Nothing, iterable) where Target + if Tables.istable(iterable) + column = Tables.getcolumn(iterable, first(GI.geometrycolumns(iterable))) + Iterators.map(x -> _flatten(f, Target, x), column) |> Iterators.flatten + else + Iterators.map(x -> _flatten(f, Target, x), iterable) |> Iterators.flatten + end +end +# Flatten feature collections +function _flatten(f, ::Type{Target}, ::GI.FeatureCollectionTrait, fc) where Target + Iterators.map(GI.getfeature(fc)) do feature + _flatten(f, Target, feature) + end |> Iterators.flatten +end +_flatten(f, ::Type{Target}, ::GI.FeatureTrait, feature) where Target = + _flatten(f, Target, GI.geometry(feature)) +# Apply f to the target geometry +_flatten(f, ::Type{Target}, ::Trait, geom) where {Target,Trait<:Target} = (f(geom),) +_flatten(f, ::Type{Target}, trait, geom) where Target = + Iterators.flatten(Iterators.map(g -> _flatten(f, Target, g), GI.getgeom(geom))) +# Fail if we hit PointTrait without running `f` +_flatten(f, ::Type{Target}, trait::GI.PointTrait, geom) where Target = + throw(ArgumentError("target $Target not found, but reached a `PointTrait` leaf")) +# Specific cases to avoid method ambiguity +_flatten(f, ::Type{<:GI.PointTrait}, ::GI.PointTrait, geom) = (f(geom),) +_flatten(f, ::Type{<:GI.FeatureTrait}, ::GI.FeatureTrait, feature) = (f(feature),) +_flatten(f, ::Type{<:GI.FeatureCollectionTrait}, ::GI.FeatureCollectionTrait, fc) = (f(fc),) + + +""" + reconstruct(geom, components) + +Reconstruct `geom` from an iterable of component objects that match its structure. + +All objects in `components` must have the same `GeoInterface.trait`. + +Usually used in combination with `flatten`. +""" +function reconstruct(geom, components) + obj, iter = _reconstruct(geom, components) + return obj +end + +_reconstruct(geom, components) = + _reconstruct(typeof(GI.trait(first(components))), geom, components, 1) +_reconstruct(::Type{Target}, geom, components, iter) where Target = + _reconstruct(Target, GI.trait(geom), geom, components, iter) +# Try to reconstruct over iterables +function _reconstruct(::Type{Target}, ::Nothing, iterable, components, iter) where Target + vect = map(iterable) do x + # iter is updated by _reconstruct here + obj, iter = _reconstruct(Target, x, components, iter) + obj + end + return vect, iter +end +# Reconstruct feature collections +function _reconstruct(::Type{Target}, ::GI.FeatureCollectionTrait, fc, components, iter) where Target + features = map(GI.getfeature(fc)) do feature + # iter is updated by _reconstruct here + newfeature, iter = _reconstruct(Target, feature, components, iter) + newfeature + end + return GI.FeatureCollection(features; crs=GI.crs(fc)), iter +end +function _reconstruct(::Type{Target}, ::GI.FeatureTrait, feature, components, iter) where Target + geom, iter = _reconstruct(Target, GI.geometry(feature), components, iter) + return GI.Feature(geom; properties=GI.properties(feature), crs=GI.crs(feature)), iter +end +function _reconstruct(::Type{Target}, trait, geom, components, iter) where Target + geoms = map(GI.getgeom(geom)) do subgeom + # iter is updated by _reconstruct here + subgeom1, iter = _reconstruct(Target, GI.trait(subgeom), subgeom, components, iter) + subgeom1 + end + return rebuild(geom, geoms), iter +end +# Apply f to the target geometry +_reconstruct(::Type{Target}, ::Trait, geom, components, iter) where {Target,Trait<:Target} = + iterate(components, iter) +# Specific cases to avoid method ambiguity +_reconstruct(::Type{<:GI.PointTrait}, ::GI.PointTrait, geom, components, iter) = iterate(components, iter) +_reconstruct(::Type{<:GI.FeatureTrait}, ::GI.FeatureTrait, feature, components, iter) = iterate(feature, iter) +_reconstruct(::Type{<:GI.FeatureCollectionTrait}, ::GI.FeatureCollectionTrait, fc, components, iter) = iterate(fc, iter) +# Fail if we hit PointTrait without running `f` +_reconstruct(::Type{Target}, trait::GI.PointTrait, geom, components, iter) where Target = + throw(ArgumentError("target $Target not found, but reached a `PointTrait` leaf")) + + +const BasicsGeoms = Union{GB.AbstractGeometry,GB.AbstractFace,GB.AbstractPoint,GB.AbstractMesh, + GB.AbstractPolygon,GB.LineString,GB.MultiPoint,GB.MultiLineString,GB.MultiPolygon,GB.Mesh} + +""" + rebuild(geom, child_geoms) + +Rebuild a geometry from child geometries. + +By default geometries will be rebuilt as a `GeoInterface.Wrappers` +geometry, but `rebuild` can have methods added to it to dispatch +on geometries from other packages and specify how to rebuild them. + +(Maybe it should go into GeoInterface.jl) +""" +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 +end diff --git a/GeometryOpsCore/src/types.jl b/GeometryOpsCore/src/types.jl new file mode 100644 index 000000000..d14a44df2 --- /dev/null +++ b/GeometryOpsCore/src/types.jl @@ -0,0 +1,95 @@ +#= +# Types + +This defines core types that the GeometryOps ecosystem uses, +and that are usable in more than just GeometryOps. + +=# + +export TraitTarget +export BoolsAsTypes, _True, _False, _booltype + +#= +## `TraitTarget` + +This struct holds a trait parameter or a union of trait parameters. +It's essentially a way to construct unions. +=# + +""" + TraitTarget{T} + +This struct holds a trait parameter or a union of trait parameters. + +It is primarily used for dispatch into methods which select trait levels, +like `apply`, or as a parameter to `target`. + +## Constructors +```julia +TraitTarget(GI.PointTrait()) +TraitTarget(GI.LineStringTrait(), GI.LinearRingTrait()) # and other traits as you may like +TraitTarget(TraitTarget(...)) +# There are also type based constructors available, but that's not advised. +TraitTarget(GI.PointTrait) +TraitTarget(Union{GI.LineStringTrait, GI.LinearRingTrait}) +# etc. +``` + +""" +struct TraitTarget{T} end +TraitTarget(::Type{T}) where T = TraitTarget{T}() +TraitTarget(::T) where T<:GI.AbstractTrait = TraitTarget{T}() +TraitTarget(::TraitTarget{T}) where T = TraitTarget{T}() +TraitTarget(::Type{<:TraitTarget{T}}) where T = TraitTarget{T}() +TraitTarget(traits::GI.AbstractTrait...) = TraitTarget{Union{map(typeof, traits)...}}() + + +Base.in(::Trait, ::TraitTarget{Target}) where {Trait <: GI.AbstractTrait, Target} = Trait <: Target + + + +#= +## `BoolsAsTypes` + +In `apply` and `applyreduce`, we pass `threading` and `calc_extent` as types, not simple boolean values. + +This is to help compilation - with a type to hold on to, it's easier for +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. + +TODO: should we switch to `Static.jl`? +=# + +""" + abstract type BoolsAsTypes + +""" +abstract type BoolsAsTypes end + +""" + struct _True <: BoolsAsTypes + +A struct that means `true`. +""" +struct _True <: BoolsAsTypes end + +""" + struct _False <: BoolsAsTypes + +A struct that means `false`. +""" +struct _False <: BoolsAsTypes end + +""" + _booltype(x) + +Returns a [`BoolsAsTypes`](@ref) from `x`, whether it's a boolean or a BoolsAsTypes. +""" +function _booltype end + +@inline _booltype(x::Bool)::BoolsAsTypes = x ? _True() : _False() +@inline _booltype(x::BoolsAsTypes)::BoolsAsTypes = x diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index 308ad7383..ff4a6b1b3 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -2,6 +2,17 @@ module GeometryOps +include("../GeometryOpsCore/src/GeometryOpsCore.jl") +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 + using GeoInterface using GeometryBasics using LinearAlgebra, Statistics diff --git a/src/primitives.jl b/src/primitives.jl index 34a315fcd..e69de29bb 100644 --- a/src/primitives.jl +++ b/src/primitives.jl @@ -1,652 +0,0 @@ -# # Primitive functions - -export apply, applyreduce, TraitTarget - -#= - -This file mainly defines the [`apply`](@ref) and [`applyreduce`](@ref) functions, and some related functionality. - -In general, the idea behind the `apply` framework is to take -as input any geometry, vector of geometries, or feature collection, -deconstruct it to the given trait target (any arbitrary GI.AbstractTrait -or `TraitTarget` union thereof, like `PointTrait` or `PolygonTrait`) -and perform some operation on it. - -This allows for a simple and consistent framework within which users can -define their own operations trivially easily, and removes a lot of the -complexity involved with handling complex geometry structures. - -For example, a simple way to flip the x and y coordinates of a geometry is: - -```julia -flipped_geom = GO.apply(GI.PointTrait(), geom) do p - (GI.y(p), GI.x(p)) -end -``` - -As simple as that. There's no need to implement your own decomposition because it's done for you. - -Functions like [`flip`](@ref), [`reproject`](@ref), [`transform`](@ref), even [`segmentize`](@ref) and [`simplify`](@ref) have been implemented -using the `apply` framework. Similarly, [`centroid`](@ref), [`area`](@ref) and [`distance`](@ref) have been implemented using the -[`applyreduce`](@ref) framework. - -## Docstrings - -### Functions - -```@docs -apply -applyreduce -GeometryOps.unwrap -GeometryOps.flatten -GeometryOps.reconstruct -GeometryOps.rebuild -``` - -## Types - -```@docs -TraitTarget -``` - -## Implementation - -=# - -const THREADED_KEYWORD = "- `threaded`: `true` or `false`. Whether to use multithreading. Defaults to `false`." -const CRS_KEYWORD = "- `crs`: The CRS to attach to geometries. Defaults to `nothing`." -const CALC_EXTENT_KEYWORD = "- `calc_extent`: `true` or `false`. Whether to calculate the extent. Defaults to `false`." - -const APPLY_KEYWORDS = """ -$THREADED_KEYWORD -$CRS_KEYWORD -$CALC_EXTENT_KEYWORD -""" - -#= -## What is `apply`? - -`apply` applies some function to every geometry matching the `Target` -GeoInterface trait, in some arbitrarily nested object made up of: -- `AbstractArray`s (we also try to iterate other non-GeoInteface compatible object) -- `FeatureCollectionTrait` objects -- `FeatureTrait` objects -- `AbstractGeometryTrait` objects - -`apply` recursively calls itself through these nested -layers until it reaches objects with the `Target` GeoInterface trait. When found `apply` applies the function `f`, and stops. - -The outer recursive functions then progressively rebuild the object -using GeoInterface objects matching the original traits. - -If `PointTrait` is found but it is not the `Target`, an error is thrown. -This likely means the object contains a different geometry trait to -the target, such as `MultiPointTrait` when `LineStringTrait` was specified. - -To handle this possibility it may be necessary to make `Target` a -`Union` of traits found at the same level of nesting, and define methods -of `f` to handle all cases. - -Be careful making a union across "levels" of nesting, e.g. -`Union{FeatureTrait,PolygonTrait}`, as `_apply` will just never reach -`PolygonTrait` when all the polygons are wrapped in a `FeatureTrait` object. - -## Embedding: - -`extent` and `crs` can be embedded in all geometries, features, and -feature collections as part of `apply`. Geometries deeper than `Target` -will of course not have new `extent` or `crs` embedded. - -- `calc_extent` signals to recalculate an `Extent` and embed it. -- `crs` will be embedded as-is - -## Threading - -Threading is used at the outermost level possible - over -an array, feature collection, or e.g. a MultiPolygonTrait where -each `PolygonTrait` sub-geometry may be calculated on a different thread. - -Currently, threading defaults to `false` for all objects, but can be turned on -by passing the keyword argument `threaded=true` to `apply`. -=# - -""" - apply(f, target::Union{TraitTarget, GI.AbstractTrait}, obj; kw...) - -Reconstruct a geometry, feature, feature collection, or nested vectors of -either using the function `f` on the `target` trait. - -`f(target_geom) => x` where `x` also has the `target` trait, or a trait that can -be substituted. For example, swapping `PolgonTrait` to `MultiPointTrait` will fail -if the outer object has `MultiPolygonTrait`, but should work if it has `FeatureTrait`. - -Objects "shallower" than the target trait are always completely rebuilt, like -a `Vector` of `FeatureCollectionTrait` of `FeatureTrait` when the target -has `PolygonTrait` and is held in the features. These will always be GeoInterface -geometries/feature/feature collections. But "deeper" objects may remain -unchanged or be whatever GeoInterface compatible objects `f` returns. - -The result is a functionally similar geometry with values depending on `f`. - -$APPLY_KEYWORDS - -## Example - -Flipped point the order in any feature or geometry, or iterables of either: - -```julia -import GeoInterface as GI -import GeometryOps as GO -geom = GI.Polygon([GI.LinearRing([(1, 2), (3, 4), (5, 6), (1, 2)]), - GI.LinearRing([(3, 4), (5, 6), (6, 7), (3, 4)])]) - -flipped_geom = GO.apply(GI.PointTrait, geom) do p - (GI.y(p), GI.x(p)) -end -``` -""" -@inline function apply( - f::F, target, geom; calc_extent=false, threaded=false, kw... -) where F - threaded = _booltype(threaded) - calc_extent = _booltype(calc_extent) - _apply(f, TraitTarget(target), geom; threaded, calc_extent, kw...) -end - -# Call _apply again with the trait of `geom` -@inline _apply(f::F, target, geom; kw...) where F = - _apply(f, target, GI.trait(geom), geom; kw...) -# There is no trait and this is an AbstractArray - so just iterate over it calling _apply on the contents -@inline function _apply(f::F, target, ::Nothing, A::AbstractArray; threaded, kw...) where F - # 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...) - _maptasks(apply_to_array, eachindex(A), threaded) -end -# There is no trait and this is not an AbstractArray. -# Try to call _apply over it. We can't use threading -# as we don't know if we can can index into it. So just `map`. -@inline function _apply(f::F, target, ::Nothing, iterable::IterableType; threaded, kw...) where {F, IterableType} - # Try the Tables.jl interface first - if Tables.istable(iterable) - _apply_table(f, target, iterable; threaded, kw...) - else # this is probably some form of iterable... - if threaded isa _True - # `collect` first so we can use threads - _apply(f, target, collect(iterable); threaded, kw...) - else - apply_to_iterable(x) = _apply(f, target, x; kw...) - map(apply_to_iterable, iterable) - end - end -end -#= -Doing this inline in `_apply` is _heavily_ type unstable, so it's best to separate this -by a function barrier. - -This function operates `apply` on the `geometry` column of the table, and returns a new table -with the same schema, but with the new geometry column. - -This new table may be of the same type as the old one iff `Tables.materializer` is defined for -that table. If not, then a `NamedTuple` is returned. -=# -function _apply_table(f::F, target, iterable::IterableType; threaded, kw...) where {F, IterableType} - _get_col_pair(colname) = colname => Tables.getcolumn(iterable, colname) - # We extract the geometry column and run `apply` on it. - geometry_column = first(GI.geometrycolumns(iterable)) - new_geometry = _apply(f, target, Tables.getcolumn(iterable, geometry_column); threaded, kw...) - # Then, we obtain the schema of the table, - old_schema = Tables.schema(iterable) - # filter the geometry column out, - new_names = filter(Base.Fix1(!==, geometry_column), old_schema.names) - # and try to rebuild the same table as the best type - either the original type of `iterable`, - # or a named tuple which is the default fallback. - result = Tables.materializer(iterable)( - merge( - NamedTuple{(geometry_column,), Base.Tuple{typeof(new_geometry)}}((new_geometry,)), - NamedTuple(Iterators.map(_get_col_pair, new_names)) - ) - ) - # Finally, we ensure that metadata is propagated correctly. - # This can only happen if the original table supports metadata reads, - # and the result supports metadata writes. - if DataAPI.metadatasupport(typeof(result)).write - # Copy over all metadata from the original table to the new table, - # if the original table supports metadata reading. - if DataAPI.metadatasupport(IterableType).read - for (key, (value, style)) in DataAPI.metadata(iterable; style = true) - # Default styles are not preserved on data transformation, so we must skip them! - style == :default && continue - # We assume that any other style is preserved. - DataAPI.metadata!(result, key, value; style) - end - end - # We don't usually care about the original table's metadata for GEOINTERFACE namespaced - # keys, so we should set the crs and geometrycolumns metadata if they are present. - # Ensure that `GEOINTERFACE:geometrycolumns` and `GEOINTERFACE:crs` are set! - mdk = DataAPI.metadatakeys(result) - # If the user has asked for geometry columns to persist, they would be here, - # so we don't need to set them. - if !("GEOINTERFACE:geometrycolumns" in mdk) - # If the geometry columns are not already set, we need to set them. - DataAPI.metadata!(result, "GEOINTERFACE:geometrycolumns", (geometry_column,); style = :default) - end - # Force reset CRS always, since you can pass `crs` to `apply`. - new_crs = if haskey(kw, :crs) - kw[:crs] - else - GI.crs(iterable) # this will automatically check `GEOINTERFACE:crs` unless the type has a specialized implementation. - end - - DataAPI.metadata!(result, "GEOINTERFACE:crs", new_crs; style = :default) - end - - return result -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 -) 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 - features = _maptasks(apply_to_feature, 1:GI.nfeature(fc), threaded) - 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 - return GI.FeatureCollection(features; crs, extent) - else - # Return a FeatureCollection with features and crs - return GI.FeatureCollection(features; crs) - 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 -) 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 - # 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 - return GI.Feature(geometry; properties, crs, extent) - else - # Return a new Feature with the new geometry, but the original properties and crs - return GI.Feature(geometry; properties, crs) - end -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 -)::(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()) - 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 -)::(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()) - 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) - # 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) - # Return a new geometry of the same trait as `geom`, holding the new `geoms` with `crs` - return rebuild(geom, geoms; crs) -end -# Fail loudly if we hit PointTrait without running `f` -# (after PointTrait there is no further to dig with `_apply`) -# @inline _apply(f, ::TraitTarget{Target}, trait::GI.PointTrait, geom; crs=nothing, kw...) where Target = - # throw(ArgumentError("target $Target not found, but reached a `PointTrait` leaf")) -# Finally, these short methods are the main purpose of `apply`. -# The `Trait` is a subtype of the `Target` (or identical to it) -# So the `Target` is found. We apply `f` to geom and return it to previous -# _apply calls to be wrapped with the outer geometries/feature/featurecollection/array. -_apply(f::F, ::TraitTarget{Target}, ::Trait, geom; crs=GI.crs(geom), kw...) where {F,Target,Trait<:Target} = f(geom) -# Define some specific cases of this match to avoid method ambiguity -for T in ( - GI.PointTrait, GI.LinearRing, GI.LineString, - GI.MultiPoint, GI.FeatureTrait, GI.FeatureCollectionTrait -) - @eval _apply(f::F, target::TraitTarget{<:$T}, trait::$T, x; kw...) where F = f(x) -end - -""" - applyreduce(f, op, target::Union{TraitTarget, GI.AbstractTrait}, obj; threaded) - -Apply function `f` to all objects with the `target` trait, -and reduce the result with an `op` like `+`. - -The order and grouping of application of `op` is not guaranteed. - -If `threaded==true` threads will be used over arrays and iterables, -feature collections and nested geometries. -""" -@inline function applyreduce( - f::F, op::O, target, geom; threaded=false, init=nothing -) where {F, O} - 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) - _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(f, op, target, collect(iterable); threaded, init) - else - # Try to `mapreduce` the iterable as-is - mapreduce(applyreduce_iterable, op, iterable; init) - end - end -end -# In this case, we don't reconstruct the table, but only operate on the geometry column. -function _applyreduce_table(f::F, op::O, target, iterable::IterableType; threaded, init) where {F, O, IterableType} - # We extract the geometry column and run `applyreduce` on it. - geometry_column = first(GI.geometrycolumns(iterable)) - return _applyreduce(f, op, target, Tables.getcolumn(iterable, geometry_column); threaded, init) -end -# If `applyreduce` wants features, then applyreduce over the rows as `GI.Feature`s. -function _applyreduce_table(f::F, op::O, target::GI.FeatureTrait, iterable::IterableType; threaded, init) where {F, O, IterableType} - # We extract the geometry column and run `apply` on it. - geometry_column = first(GI.geometrycolumns(iterable)) - property_names = Iterators.filter(!=(geometry_column), Tables.schema(iterable).names) - features = map(Tables.rows(iterable)) do row - GI.Feature(Tables.getcolumn(row, geometry_column), properties=NamedTuple(Iterators.map(Base.Fix1(_get_col_pair, row), property_names))) - end - return _applyreduce(f, op, target, features; threaded, init) -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) - _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) - _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) -end -# Apply f to the target -@inline function _applyreduce(f::F, op::O, ::TraitTarget{Target}, ::Trait, x; kw...) where {F,O,Target,Trait<:Target} - f(x) -end -# Fail if we hit PointTrait -# _applyreduce(f, op, target::TraitTarget{Target}, trait::PointTrait, geom; kw...) where Target = - # throw(ArgumentError("target $target not found")) -# Specific cases to avoid method ambiguity -for T in ( - GI.PointTrait, GI.LinearRing, GI.LineString, - GI.MultiPoint, GI.FeatureTrait, GI.FeatureCollectionTrait -) - @eval _applyreduce(f::F, op::O, ::TraitTarget{<:$T}, trait::$T, x; kw...) where {F, O} = f(x) -end - -""" - unwrap(target::Type{<:AbstractTrait}, obj) - unwrap(f, target::Type{<:AbstractTrait}, obj) - -Unwrap the object to vectors, down to the target trait. - -If `f` is passed in it will be applied to the target geometries -as they are found. -""" -function unwrap end -unwrap(target::Type, geom) = unwrap(identity, target, geom) -# Add dispatch argument for trait -unwrap(f, target::Type, geom) = unwrap(f, target, GI.trait(geom), geom) -# Try to unwrap over iterables -unwrap(f, target::Type, ::Nothing, iterable) = - map(x -> unwrap(f, target, x), iterable) -# Rewrap feature collections -unwrap(f, target::Type, ::GI.FeatureCollectionTrait, fc) = - map(x -> unwrap(f, target, x), GI.getfeature(fc)) -unwrap(f, target::Type, ::GI.FeatureTrait, feature) = - unwrap(f, target, GI.geometry(feature)) -unwrap(f, target::Type, trait, geom) = map(g -> unwrap(f, target, g), GI.getgeom(geom)) -# Apply f to the target geometry -unwrap(f, ::Type{Target}, ::Trait, geom) where {Target,Trait<:Target} = f(geom) -# Fail if we hit PointTrait -unwrap(f, target::Type, trait::GI.PointTrait, geom) = - throw(ArgumentError("target $target not found, but reached a `PointTrait` leaf")) -# Specific cases to avoid method ambiguity -unwrap(f, target::Type{GI.PointTrait}, trait::GI.PointTrait, geom) = f(geom) -unwrap(f, target::Type{GI.FeatureTrait}, ::GI.FeatureTrait, feature) = f(feature) -unwrap(f, target::Type{GI.FeatureCollectionTrait}, ::GI.FeatureCollectionTrait, fc) = f(fc) - -""" - flatten(target::Type{<:GI.AbstractTrait}, obj) - flatten(f, target::Type{<:GI.AbstractTrait}, obj) - -Lazily flatten any `AbstractArray`, iterator, `FeatureCollectionTrait`, -`FeatureTrait` or `AbstractGeometryTrait` object `obj`, so that objects -with the `target` trait are returned by the iterator. - -If `f` is passed in it will be applied to the target geometries. -""" -flatten(::Type{Target}, geom) where {Target<:GI.AbstractTrait} = flatten(identity, Target, geom) -flatten(f, ::Type{Target}, geom) where {Target<:GI.AbstractTrait} = _flatten(f, Target, geom) - -_flatten(f, ::Type{Target}, geom) where Target = _flatten(f, Target, GI.trait(geom), geom) -# Try to flatten over iterables -function _flatten(f, ::Type{Target}, ::Nothing, iterable) where Target - if Tables.istable(iterable) - column = Tables.getcolumn(iterable, first(GI.geometrycolumns(iterable))) - Iterators.map(x -> _flatten(f, Target, x), column) |> Iterators.flatten - else - Iterators.map(x -> _flatten(f, Target, x), iterable) |> Iterators.flatten - end -end -# Flatten feature collections -function _flatten(f, ::Type{Target}, ::GI.FeatureCollectionTrait, fc) where Target - Iterators.map(GI.getfeature(fc)) do feature - _flatten(f, Target, feature) - end |> Iterators.flatten -end -_flatten(f, ::Type{Target}, ::GI.FeatureTrait, feature) where Target = - _flatten(f, Target, GI.geometry(feature)) -# Apply f to the target geometry -_flatten(f, ::Type{Target}, ::Trait, geom) where {Target,Trait<:Target} = (f(geom),) -_flatten(f, ::Type{Target}, trait, geom) where Target = - Iterators.flatten(Iterators.map(g -> _flatten(f, Target, g), GI.getgeom(geom))) -# Fail if we hit PointTrait without running `f` -_flatten(f, ::Type{Target}, trait::GI.PointTrait, geom) where Target = - throw(ArgumentError("target $Target not found, but reached a `PointTrait` leaf")) -# Specific cases to avoid method ambiguity -_flatten(f, ::Type{<:GI.PointTrait}, ::GI.PointTrait, geom) = (f(geom),) -_flatten(f, ::Type{<:GI.FeatureTrait}, ::GI.FeatureTrait, feature) = (f(feature),) -_flatten(f, ::Type{<:GI.FeatureCollectionTrait}, ::GI.FeatureCollectionTrait, fc) = (f(fc),) - - -""" - reconstruct(geom, components) - -Reconstruct `geom` from an iterable of component objects that match its structure. - -All objects in `components` must have the same `GeoInterface.trait`. - -Usually used in combination with `flatten`. -""" -function reconstruct(geom, components) - obj, iter = _reconstruct(geom, components) - return obj -end - -_reconstruct(geom, components) = - _reconstruct(typeof(GI.trait(first(components))), geom, components, 1) -_reconstruct(::Type{Target}, geom, components, iter) where Target = - _reconstruct(Target, GI.trait(geom), geom, components, iter) -# Try to reconstruct over iterables -function _reconstruct(::Type{Target}, ::Nothing, iterable, components, iter) where Target - vect = map(iterable) do x - # iter is updated by _reconstruct here - obj, iter = _reconstruct(Target, x, components, iter) - obj - end - return vect, iter -end -# Reconstruct feature collections -function _reconstruct(::Type{Target}, ::GI.FeatureCollectionTrait, fc, components, iter) where Target - features = map(GI.getfeature(fc)) do feature - # iter is updated by _reconstruct here - newfeature, iter = _reconstruct(Target, feature, components, iter) - newfeature - end - return GI.FeatureCollection(features; crs=GI.crs(fc)), iter -end -function _reconstruct(::Type{Target}, ::GI.FeatureTrait, feature, components, iter) where Target - geom, iter = _reconstruct(Target, GI.geometry(feature), components, iter) - return GI.Feature(geom; properties=GI.properties(feature), crs=GI.crs(feature)), iter -end -function _reconstruct(::Type{Target}, trait, geom, components, iter) where Target - geoms = map(GI.getgeom(geom)) do subgeom - # iter is updated by _reconstruct here - subgeom1, iter = _reconstruct(Target, GI.trait(subgeom), subgeom, components, iter) - subgeom1 - end - return rebuild(geom, geoms), iter -end -# Apply f to the target geometry -_reconstruct(::Type{Target}, ::Trait, geom, components, iter) where {Target,Trait<:Target} = - iterate(components, iter) -# Specific cases to avoid method ambiguity -_reconstruct(::Type{<:GI.PointTrait}, ::GI.PointTrait, geom, components, iter) = iterate(components, iter) -_reconstruct(::Type{<:GI.FeatureTrait}, ::GI.FeatureTrait, feature, components, iter) = iterate(feature, iter) -_reconstruct(::Type{<:GI.FeatureCollectionTrait}, ::GI.FeatureCollectionTrait, fc, components, iter) = iterate(fc, iter) -# Fail if we hit PointTrait without running `f` -_reconstruct(::Type{Target}, trait::GI.PointTrait, geom, components, iter) where Target = - throw(ArgumentError("target $Target not found, but reached a `PointTrait` leaf")) - - -const BasicsGeoms = Union{GB.AbstractGeometry,GB.AbstractFace,GB.AbstractPoint,GB.AbstractMesh, - GB.AbstractPolygon,GB.LineString,GB.MultiPoint,GB.MultiLineString,GB.MultiPolygon,GB.Mesh} - -""" - rebuild(geom, child_geoms) - -Rebuild a geometry from child geometries. - -By default geometries will be rebuilt as a `GeoInterface.Wrappers` -geometry, but `rebuild` can have methods added to it to dispatch -on geometries from other packages and specify how to rebuild them. - -(Maybe it should go into GeoInterface.jl) -""" -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 -end - -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 - ntasks = length(taskrange) - # Customize this as needed. - # More tasks have more overhead, but better load balancing - tasks_per_thread = 2 - chunk_size = max(1, ntasks ÷ (tasks_per_thread * nthreads())) - # partition the range into chunks - task_chunks = Iterators.partition(taskrange, chunk_size) - # Map over the chunks - tasks = map(task_chunks) do chunk - # Spawn a task to process this chunk - @spawn begin - # Where we map `f` over the chunk indices - map(f, chunk) - end - end - - # Finally we join the results into a new vector - return mapreduce(fetch, vcat, tasks) -end -#= -Here we 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! -=# -Base.@assume_effects :foldable @inline function _maptasks(f::F, taskrange, threaded::_False)::Vector where F - map(f, taskrange) -end - -# Threading utility, modified Mason Protters threading PSA -# run `f` over ntasks, where f receives an AbstractArray/range -# of linear indices -# -# WARNING: this will not work for mean/median - only ops -# where grouping is possible -@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 - tasks_per_thread = 2 - chunk_size = max(1, ntasks ÷ (tasks_per_thread * nthreads())) - # partition the range into chunks - task_chunks = Iterators.partition(taskrange, chunk_size) - # Map over the chunks - tasks = map(task_chunks) do chunk - # Spawn a task to process this chunk - @spawn begin - # Where we map `f` over the chunk indices - mapreduce(f, op, chunk; init) - end - end - - # 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 - mapreduce(f, op, taskrange; init) -end diff --git a/src/types.jl b/src/types.jl index 4e6a0dec1..3a89f0f9f 100644 --- a/src/types.jl +++ b/src/types.jl @@ -8,66 +8,7 @@ This file defines some fundamental types used in GeometryOps. This is because we define types in the files where they are used, to make it easier to understand the code. =# -export TraitTarget, GEOS -#= -## `TraitTarget` - -This struct holds a trait parameter or a union of trait parameters. -It's essentially a way to construct unions. -=# - -""" - TraitTarget{T} - -This struct holds a trait parameter or a union of trait parameters. - -It is primarily used for dispatch into methods which select trait levels, -like `apply`, or as a parameter to `target`. - -## Constructors -```julia -TraitTarget(GI.PointTrait()) -TraitTarget(GI.LineStringTrait(), GI.LinearRingTrait()) # and other traits as you may like -TraitTarget(TraitTarget(...)) -# There are also type based constructors available, but that's not advised. -TraitTarget(GI.PointTrait) -TraitTarget(Union{GI.LineStringTrait, GI.LinearRingTrait}) -# etc. -``` - -""" -struct TraitTarget{T} end -TraitTarget(::Type{T}) where T = TraitTarget{T}() -TraitTarget(::T) where T<:GI.AbstractTrait = TraitTarget{T}() -TraitTarget(::TraitTarget{T}) where T = TraitTarget{T}() -TraitTarget(::Type{<:TraitTarget{T}}) where T = TraitTarget{T}() -TraitTarget(traits::GI.AbstractTrait...) = TraitTarget{Union{map(typeof, traits)...}}() - - -Base.in(::Trait, ::TraitTarget{Target}) where {Trait <: GI.AbstractTrait, Target} = Trait <: Target - -#= -## `BoolsAsTypes` - -In `apply` and `applyreduce`, we pass `threading` and `calc_extent` as types, not simple boolean values. - -This is to help compilation - with a type to hold on to, it's easier for -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. - -TODO: should we switch to `Static.jl`? -=# -abstract type BoolsAsTypes end -struct _True <: BoolsAsTypes end -struct _False <: BoolsAsTypes end - -@inline _booltype(x::Bool)::BoolsAsTypes = x ? _True() : _False() -@inline _booltype(x::BoolsAsTypes)::BoolsAsTypes = x - +export GEOS #= ## `GEOS` From ba338a11662c587a06da2e33d10cc6e6f4afbfe7 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 2 Oct 2024 15:20:03 -0700 Subject: [PATCH 02/10] Define `Manifold` and subtypes Linear, Spherical, Geodesic The name `Manifold` is up for debate, I just wanted some supertype for dispatch constraints. --- GeometryOpsCore/src/types.jl | 64 ++++++++++++++++++++++++++++++++++++ 1 file changed, 64 insertions(+) diff --git a/GeometryOpsCore/src/types.jl b/GeometryOpsCore/src/types.jl index d14a44df2..b8b0fde88 100644 --- a/GeometryOpsCore/src/types.jl +++ b/GeometryOpsCore/src/types.jl @@ -6,9 +6,73 @@ and that are usable in more than just GeometryOps. =# +#= +## `Manifold` + +A manifold is mathematically defined as a topological space that resembles Euclidean space locally. + +In GeometryOps (and geodesy more generally), there are three manifolds we care about: +- [`Linear`](@ref): the 2d plane, a completely Euclidean manifold +- [`Spherical`](@ref): the unit sphere, but one where areas are multiplied by the radius of the Earth. This is not Euclidean globally, but all map projections attempt to represent the sphere on the Euclidean 2D plane to varying degrees of success. +- [`Geodesic`](@ref): the ellipsoid, the closest we can come to representing the Earth by a simple geometric shape. Parametrized by `semimajor_axis` and `inv_flattening`. + +Generally, we aim to have `Linear` and `Spherical` be operable everywhere, whereas `Geodesic` will only apply in specific circumstances. +Currently, those circumstances are `area` and `segmentize`, but this could be extended with time and https://github.com/JuliaGeo/SphericalGeodesics.jl. +=# + +export Linear, Spherical, Geodesic export TraitTarget export BoolsAsTypes, _True, _False, _booltype +""" + abstract type Manifold + +A manifold is mathematically defined as a topological space that resembles Euclidean space locally. + +We use the manifold definition to define the space in which an operation should be performed, or where a geometry lies. + +Currently we have [`Linear`](@ref), [`Spherical`](@ref), and [`Geodesic`](@ref) manifolds. +""" +abstract type Manifold end + +""" + Linear() + +A linear manifold means that the space is completely Euclidean, +and planar geometry suffices. +""" +struct Linear <: Manifold +end + +""" + Spherical(; radius) + +A spherical manifold means that the geometry is on the 3-sphere (but is represented by 2-D longitude and latitude). + +## Extended help + +!!! note + The traditional definition of spherical coordinates in physics and mathematics, ``r, \theta, \phi``, uses the _colatitude_, + that measures angular displacement from the `z`-axis. Here, we use the geographic definition of longitude and latitude, meaning + that `lon` is longitude between -180 and 180, and `lat` is latitude between `-90` (south pole) and `90` (north pole). +""" +Base.@kwdef struct Spherical{T} <: Manifold + radius::T = 6378137.0 +end + +""" + Geodesic(; semimajor_axis, inv_flattening) + +A geodesic manifold means that the geometry is on a 3-dimensional ellipsoid, parameterized by `semimajor_axis` (``a`` in mathematical parlance) +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 + inv_flattening::T = 298.257223563 +end + #= ## `TraitTarget` From b279306227e8935d9701a5c3e8376f5269565e81 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 2 Oct 2024 15:54:54 -0700 Subject: [PATCH 03/10] Fix integration issues, get things working --- GeometryOpsCore/src/GeometryOpsCore.jl | 11 +++++++++++ GeometryOpsCore/src/geometry_utils.jl | 3 +++ GeometryOpsCore/src/other_primitives.jl | 4 ---- GeometryOpsCore/src/types.jl | 9 ++++++--- ext/GeometryOpsLibGEOSExt/GeometryOpsLibGEOSExt.jl | 6 +++--- src/GeometryOps.jl | 6 +++--- src/utils.jl | 3 --- 7 files changed, 26 insertions(+), 16 deletions(-) create mode 100644 GeometryOpsCore/src/geometry_utils.jl diff --git a/GeometryOpsCore/src/GeometryOpsCore.jl b/GeometryOpsCore/src/GeometryOpsCore.jl index 4084f13db..bdd8fdf76 100644 --- a/GeometryOpsCore/src/GeometryOpsCore.jl +++ b/GeometryOpsCore/src/GeometryOpsCore.jl @@ -2,7 +2,17 @@ module GeometryOpsCore using Base.Threads: nthreads, @threads, @spawn +import GeoInterface import GeoInterface as GI +import GeoInterface: Extents + +# Import all names from GeoInterface and Extents, so users can do `GO.extent` or `GO.trait`. +for name in names(GeoInterface) + @eval using GeoInterface: $name +end +for name in names(Extents) + @eval using GeoInterface.Extents: $name +end using Tables using DataAPI @@ -13,5 +23,6 @@ include("types.jl") include("apply.jl") include("applyreduce.jl") include("other_primitives.jl") +include("geometry_utils.jl") end \ No newline at end of file diff --git a/GeometryOpsCore/src/geometry_utils.jl b/GeometryOpsCore/src/geometry_utils.jl new file mode 100644 index 000000000..c941273b2 --- /dev/null +++ b/GeometryOpsCore/src/geometry_utils.jl @@ -0,0 +1,3 @@ + +_linearring(geom::GI.LineString) = GI.LinearRing(parent(geom); extent=geom.extent, crs=geom.crs) +_linearring(geom::GI.LinearRing) = geom diff --git a/GeometryOpsCore/src/other_primitives.jl b/GeometryOpsCore/src/other_primitives.jl index b0fd04de7..b0ac8c6b7 100644 --- a/GeometryOpsCore/src/other_primitives.jl +++ b/GeometryOpsCore/src/other_primitives.jl @@ -145,10 +145,6 @@ _reconstruct(::Type{<:GI.FeatureCollectionTrait}, ::GI.FeatureCollectionTrait, f _reconstruct(::Type{Target}, trait::GI.PointTrait, geom, components, iter) where Target = throw(ArgumentError("target $Target not found, but reached a `PointTrait` leaf")) - -const BasicsGeoms = Union{GB.AbstractGeometry,GB.AbstractFace,GB.AbstractPoint,GB.AbstractMesh, - GB.AbstractPolygon,GB.LineString,GB.MultiPoint,GB.MultiLineString,GB.MultiPolygon,GB.Mesh} - """ rebuild(geom, child_geoms) diff --git a/GeometryOpsCore/src/types.jl b/GeometryOpsCore/src/types.jl index b8b0fde88..c43f6a97c 100644 --- a/GeometryOpsCore/src/types.jl +++ b/GeometryOpsCore/src/types.jl @@ -52,9 +52,12 @@ A spherical manifold means that the geometry is on the 3-sphere (but is represen ## Extended help !!! note - The traditional definition of spherical coordinates in physics and mathematics, ``r, \theta, \phi``, uses the _colatitude_, - that measures angular displacement from the `z`-axis. Here, we use the geographic definition of longitude and latitude, meaning - that `lon` is longitude between -180 and 180, and `lat` is latitude between `-90` (south pole) and `90` (north pole). + The traditional definition of spherical coordinates in physics and mathematics, + ``r, \\theta, \\phi``, uses the _colatitude_, that measures angular displacement from the `z`-axis. + + Here, we use the geographic definition of longitude and latitude, meaning + that `lon` is longitude between -180 and 180, and `lat` is latitude between + `-90` (south pole) and `90` (north pole). """ Base.@kwdef struct Spherical{T} <: Manifold radius::T = 6378137.0 diff --git a/ext/GeometryOpsLibGEOSExt/GeometryOpsLibGEOSExt.jl b/ext/GeometryOpsLibGEOSExt/GeometryOpsLibGEOSExt.jl index f47b8ee9b..c8cbfc9f8 100644 --- a/ext/GeometryOpsLibGEOSExt/GeometryOpsLibGEOSExt.jl +++ b/ext/GeometryOpsLibGEOSExt/GeometryOpsLibGEOSExt.jl @@ -1,7 +1,7 @@ module GeometryOpsLibGEOSExt import GeometryOps as GO, LibGEOS as LG -import GeometryOps: GI +import GeoInterface as GI import GeometryOps: GEOS, enforce @@ -11,8 +11,8 @@ using GeometryOps # However, if you import those from another module (which you would with `all=true`), # that creates an ambiguity which causes a warning during precompile/load time. # In order to avoid this, we filter out these special functions. -for name in filter(!in((:var"#eval", :eval, :var"#include", :include)), names(GeometryOps; all = true)) - @eval using GeometryOps: $name +for name in filter(!in((:var"#eval", :eval, :var"#include", :include)), names(GeometryOps)) + @eval import GeometryOps: $name end """ diff --git a/src/GeometryOps.jl b/src/GeometryOps.jl index ff4a6b1b3..f524fd87d 100644 --- a/src/GeometryOps.jl +++ b/src/GeometryOps.jl @@ -2,14 +2,14 @@ module GeometryOps -include("../GeometryOpsCore/src/GeometryOpsCore.jl") +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 + @eval import .GeometryOpsCore: $name # Re-export all exported symbols if Base.isexported(GeometryOpsCore, name) - @eval export name + @eval export $name end end diff --git a/src/utils.jl b/src/utils.jl index bf815db2e..82a1e664d 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -126,6 +126,3 @@ function _point_in_extent(p, extent::Extents.Extent) (x1, x2), (y1, y2) = extent.X, extent.Y return x1 ≤ GI.x(p) ≤ x2 && y1 ≤ GI.y(p) ≤ y2 end - -_linearring(geom::GI.LineString) = GI.LinearRing(parent(geom); extent=geom.extent, crs=geom.crs) -_linearring(geom::GI.LinearRing) = geom From 13bc96bfa5d545d6818d212994d15278cc45f0d1 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 2 Oct 2024 16:05:13 -0700 Subject: [PATCH 04/10] Fix import error --- GeometryOpsCore/src/applyreduce.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GeometryOpsCore/src/applyreduce.jl b/GeometryOpsCore/src/applyreduce.jl index 5e21c3a2d..62d4dfaa8 100644 --- a/GeometryOpsCore/src/applyreduce.jl +++ b/GeometryOpsCore/src/applyreduce.jl @@ -112,7 +112,7 @@ end ### `_mapreducetasks` - flexible, threaded mapreduce -import Base: nthreads, @threads, @spawn +import Base.Threads: nthreads, @threads, @spawn # Threading utility, modified Mason Protters threading PSA # run `f` over ntasks, where f receives an AbstractArray/range From c7a41bdefa54f865fea2c3cb0a68f55534fc0c46 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 2 Oct 2024 16:05:20 -0700 Subject: [PATCH 05/10] Finish Project.toml and README --- GeometryOpsCore/Project.toml | 5 +++++ GeometryOpsCore/README.md | 8 ++++++++ 2 files changed, 13 insertions(+) diff --git a/GeometryOpsCore/Project.toml b/GeometryOpsCore/Project.toml index 2b9dcca14..459ac2c83 100644 --- a/GeometryOpsCore/Project.toml +++ b/GeometryOpsCore/Project.toml @@ -1,3 +1,8 @@ +name = "GeometryOpsCore" +uuid = "05efe853-fabf-41c8-927e-7063c8b9f013" +authors = ["Anshul Singhvi and contributors"] +version = "0.1.0" + [deps] DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f" diff --git a/GeometryOpsCore/README.md b/GeometryOpsCore/README.md index e69de29bb..adf5a85a6 100644 --- a/GeometryOpsCore/README.md +++ b/GeometryOpsCore/README.md @@ -0,0 +1,8 @@ +# GeometryOpsCore + +This is a "core" package for [GeometryOps.jl](https://github.com/JuliaGeo/GeometryOps.jl), that defines some basic primitive functions and types for GeometryOps. + +Generally, you would depend on this to use either the GeometryOps types (like `Linear`, `Spherical`, etc) or the primitive functions like `apply`, `applyreduce`, `flatten`, etc. +All of these are also accessible from GeometryOps, so it's preferable that you use GeometryOps directly. + +Tests are in the main GeometryOps tests, we don't have separate tests for GeometryOpsCore since it's in a monorepo structure. \ No newline at end of file From 57497607edd7d80d28e054ce976fab85d4595ff1 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 2 Oct 2024 16:07:22 -0700 Subject: [PATCH 06/10] Get docs to cover core + extensions also --- docs/make.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index 80f6b6246..d10ee03a5 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -69,6 +69,8 @@ withenv("JULIA_DEBUG" => "Literate") do # allow Literate debug output to escape global literate_pages vec = [] process_literate_recursive!(vec, source_path) + process_literate_recursive!(vec, joinpath(dirname(@__DIR__), "GeometryOpsCore")) + process_literate_recursive!(vec, joinpath(dirname(@__DIR__), "ext"))s literate_pages = vec[1][2] # this is a hack to get the pages in the correct order, without an initial "src" folder. # TODO: We should probably fix the above in `process_literate_recursive!`. end @@ -78,7 +80,7 @@ download("https://hackmd.io/kpIqAR8YRJOZQDJjUKVAUQ/download", joinpath(@__DIR__, # Finally, make the docs! makedocs(; - modules=[GeometryOps], + modules=[GeometryOps, GeometryOps.GeometryOpsCore], authors="Anshul Singhvi and contributors", repo="https://github.com/JuliaGeo/GeometryOps.jl/blob/{commit}{path}#{line}", sitename="GeometryOps.jl", From f4ef996d0727ffa406d3170d1763e30328cc23a7 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 2 Oct 2024 16:13:11 -0700 Subject: [PATCH 07/10] Fix typo --- docs/make.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index d10ee03a5..e8c530050 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -70,7 +70,7 @@ withenv("JULIA_DEBUG" => "Literate") do # allow Literate debug output to escape vec = [] process_literate_recursive!(vec, source_path) process_literate_recursive!(vec, joinpath(dirname(@__DIR__), "GeometryOpsCore")) - process_literate_recursive!(vec, joinpath(dirname(@__DIR__), "ext"))s + process_literate_recursive!(vec, joinpath(dirname(@__DIR__), "ext")) literate_pages = vec[1][2] # this is a hack to get the pages in the correct order, without an initial "src" folder. # TODO: We should probably fix the above in `process_literate_recursive!`. end From 4e6fa215591985ee236bbe8871276b1bea416367 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 2 Oct 2024 16:27:37 -0700 Subject: [PATCH 08/10] push things properly --- docs/make.jl | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index e8c530050..c4c9aa096 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -69,9 +69,15 @@ withenv("JULIA_DEBUG" => "Literate") do # allow Literate debug output to escape global literate_pages vec = [] process_literate_recursive!(vec, source_path) - process_literate_recursive!(vec, joinpath(dirname(@__DIR__), "GeometryOpsCore")) - process_literate_recursive!(vec, joinpath(dirname(@__DIR__), "ext")) literate_pages = vec[1][2] # this is a hack to get the pages in the correct order, without an initial "src" folder. + + core_stuff = [] + process_literate_recursive!(core_stuff, joinpath(dirname(@__DIR__), "GeometryOpsCore")) + push!(literate_pages, "GeometryOpsCore" => core_stuff[1][2]) + + ext_stuff = [] + process_literate_recursive!(ext_stuff, joinpath(dirname(@__DIR__), "ext")) + push!(literate_pages, "Extensions" => ext_stuff[1][2]) # TODO: We should probably fix the above in `process_literate_recursive!`. end From 35972f0b859d9f5b308bb5f2bfbba0b630a369c7 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 2 Oct 2024 17:19:00 -0700 Subject: [PATCH 09/10] amend paths too --- docs/make.jl | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index c4c9aa096..c8e6091fa 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -70,11 +70,15 @@ withenv("JULIA_DEBUG" => "Literate") do # allow Literate debug output to escape vec = [] process_literate_recursive!(vec, source_path) literate_pages = vec[1][2] # this is a hack to get the pages in the correct order, without an initial "src" folder. - + + global source_path + source_path = joinpath(dirname(@__DIR__), "GeometryOpsCore") core_stuff = [] - process_literate_recursive!(core_stuff, joinpath(dirname(@__DIR__), "GeometryOpsCore")) + process_literate_recursive!(core_stuff, ) push!(literate_pages, "GeometryOpsCore" => core_stuff[1][2]) + global source_path + source_path = joinpath(dirname(@__DIR__), "ext") ext_stuff = [] process_literate_recursive!(ext_stuff, joinpath(dirname(@__DIR__), "ext")) push!(literate_pages, "Extensions" => ext_stuff[1][2]) From 2fa161cea5a1c3621352e274c515cc2f4a2dd439 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Wed, 2 Oct 2024 17:24:46 -0700 Subject: [PATCH 10/10] fix again --- docs/make.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index c8e6091fa..fb2ac0baf 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -74,13 +74,13 @@ withenv("JULIA_DEBUG" => "Literate") do # allow Literate debug output to escape global source_path source_path = joinpath(dirname(@__DIR__), "GeometryOpsCore") core_stuff = [] - process_literate_recursive!(core_stuff, ) + process_literate_recursive!(core_stuff, source_path) push!(literate_pages, "GeometryOpsCore" => core_stuff[1][2]) global source_path source_path = joinpath(dirname(@__DIR__), "ext") ext_stuff = [] - process_literate_recursive!(ext_stuff, joinpath(dirname(@__DIR__), "ext")) + process_literate_recursive!(ext_stuff, source_path) push!(literate_pages, "Extensions" => ext_stuff[1][2]) # TODO: We should probably fix the above in `process_literate_recursive!`. end