Merge branch 'main' into as/trees

Author: asinghvi17

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.5"
+version = "0.1.6"
 
 [deps]
 DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"

GeometryOpsCore/src/apply.jl

@@ -175,9 +175,11 @@ 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; geometrycolumn = nothing, preserve_default_metadata = false, threaded, kw...) where {F, IterableType}
-    _get_col_pair(colname) = colname => Tables.getcolumn(iterable, colname)
-    # We extract the geometry column and run `apply` on it.
+function _apply_table(f::F, target, iterable::IterableType; geometrycolumn = nothing, preserve_default_metadata = false, threaded, kw...) where {F, IterableType}    # We extract the geometry column and run `apply` on it.
+    # First, we need the table schema:
+    input_schema = Tables.schema(iterable)
+    input_colnames = input_schema.names
+    # then, we find the geometry column(s)
     geometry_columns = if isnothing(geometrycolumn)
         GI.geometrycolumns(iterable)
     elseif geometrycolumn isa NTuple{N, <: Symbol} where N
@@ -187,31 +189,31 @@ function _apply_table(f::F, target, iterable::IterableType; geometrycolumn = not
     else
         throw(ArgumentError("geometrycolumn must be a Symbol or a tuple of Symbols, got a $(typeof(geometrycolumn))"))
     end
-    if !all(Base.Fix2(in, Tables.columnnames(iterable)), geometry_columns)
+    if !Base.issubset(geometry_columns, input_colnames)
         throw(ArgumentError(
             """
             `apply`: the `geometrycolumn` kwarg must be a subset of the column names of the table, 
             got $(geometry_columns)
             but the table has columns 
-            $(Tables.columnnames(iterable))
+            $(input_colnames)
             """
             ))
     end
+    # here we apply the function to the geometry column(s).
+    apply_kw = if isempty(used_reconstruct_table_kwargs(iterable))
+        kw
+    else
+        Base.structdiff(values(kw), NamedTuple{used_reconstruct_table_kwargs(iterable)})
+    end
     new_geometry_vecs = map(geometry_columns) do colname
-        _apply(f, target, Tables.getcolumn(iterable, colname); threaded, kw...)
+        _apply(f, target, Tables.getcolumn(iterable, colname); threaded, apply_kw...)
     end
-    # Then, we obtain the schema of the table,
-    old_schema = Tables.schema(iterable)
-    # filter the geometry column out,
-    new_names = filter(x -> !(x in geometry_columns), old_schema.names)
+    # Then, we filter the geometry column(s) out,
+    new_names = filter(x -> !(x in geometry_columns), input_colnames)
     # 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_columns, Base.Tuple{typeof.(new_geometry_vecs)...}}(new_geometry_vecs),
-            NamedTuple(Iterators.map(_get_col_pair, new_names))
-        )
-    )
+    # See the function directly below this one for the actual fallback implementation.
+    result = reconstruct_table(iterable, geometry_columns, new_geometry_vecs, new_names; kw...)
     # 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.
@@ -246,6 +248,51 @@ function _apply_table(f::F, target, iterable::IterableType; geometrycolumn = not
     return result
 end
 
+
+"""
+    used_reconstruct_table_kwargs(input)
+
+Return a tuple of the kwargs that should be passed to `reconstruct_table` for the given input.
+
+This is "semi-public" API, and required for any input type that defines `reconstruct_table`.
+"""
+function used_reconstruct_table_kwargs(input) 
+    ()
+end
+
+"""
+    reconstruct_table(input, geometry_column_names, geometry_columns, other_column_names, args...; kwargs...)
+
+Reconstruct a table from the given input, geometry column names,
+geometry columns, and other column names.
+
+Any function that defines `reconstruct_table` must also define `used_reconstruct_table_kwargs`.
+
+The input must be a table.
+
+The function should return a best-effort attempt at a table of the same type as the input,
+with the new geometry column(s) and other columns.
+
+The fallback implementation invokes `Tables.materializer`.  But if you want to be efficient 
+and pass e.g. arbitrary kwargs to the materializer, or materialize in a different way, you 
+can do so by overloading this function for your desired input type.
+
+This is "semi-public" API and while it may add optional arguments, it will not add new required 
+positional arguments. All implementations must allow arbitrary kwargs to pass through and harvest 
+what they need.
+"""
+function reconstruct_table(input, geometry_column_names, geometry_columns, other_column_names, args...; kwargs...)
+    @assert Tables.istable(input)
+    _get_col_pair(colname) = colname => Tables.getcolumn(input, colname)
+
+    return Tables.materializer(input)(
+        merge(
+            NamedTuple{geometry_column_names, Base.Tuple{typeof.(geometry_columns)...}}(geometry_columns),
+            NamedTuple(Iterators.map(_get_col_pair, other_column_names))
+        )
+    )
+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;
@@ -323,12 +370,20 @@ end
 # 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)
+function _apply(a::WithTrait{F}, ::TraitTarget{Target}, trait::Trait, geom; crs=GI.crs(geom), kw...) where {F,Target,Trait<:Target} 
+    a(trait, geom; Base.structdiff(values(kw), NamedTuple{(:threaded, :calc_extent)})...)
+end
 # 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)
+    @eval begin 
+        _apply(f::F, target::TraitTarget{<:$T}, trait::$T, x; kw...) where F = f(x)
+        function _apply(a::WithTrait{F}, target::TraitTarget{<:$T}, trait::$T, x; kw...) where F
+            a(trait, x; Base.structdiff(values(kw), NamedTuple{(:threaded, :calc_extent)})...)
+        end
+    end
 end
 
 

GeometryOpsCore/src/applyreduce.jl

@@ -51,7 +51,7 @@ Literate.jl source code is below.
 =#
 
 """
-    applyreduce(f, op, target::Union{TraitTarget, GI.AbstractTrait}, obj; threaded)
+    applyreduce(f, op, target::Union{TraitTarget, GI.AbstractTrait}, obj; threaded, init, kw...)
 
 Apply function `f` to all objects with the `target` trait,
 and reduce the result with an `op` like `+`. 
@@ -60,6 +60,8 @@ 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.
+
+`init` functions the same way as it does in base Julia functions like `reduce`.
 """
 @inline function applyreduce(
     f::F, op::O, target, geom; threaded=false, init=nothing
@@ -129,6 +131,9 @@ end
 @inline function _applyreduce(f::F, op::O, ::TraitTarget{Target}, ::Trait, x; kw...) where {F,O,Target,Trait<:Target} 
     f(x)
 end
+@inline function _applyreduce(a::WithTrait{F}, op::O, ::TraitTarget{Target}, trait::Trait, x; kw...) where {F,O,Target,Trait<:Target} 
+    a(trait, x; Base.structdiff(values(kw), NamedTuple{(:threaded, :init)})...)
+end
 # Fail if we hit PointTrait
 # _applyreduce(f, op, target::TraitTarget{Target}, trait::PointTrait, geom; kw...) where Target = 
     # throw(ArgumentError("target $target not found"))
@@ -137,7 +142,12 @@ 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)
+    @eval begin
+        _applyreduce(f::F, op::O, ::TraitTarget{<:$T}, trait::$T, x; kw...) where {F, O} = f(x)
+        function _applyreduce(a::WithTrait{F}, op::O, ::TraitTarget{<:$T}, trait::$T, x; kw...) where {F, O}
+            a(trait, x; Base.structdiff(values(kw), NamedTuple{(:threaded, :init)})...)
+        end
+    end
 end
 
 ### `_mapreducetasks` - flexible, threaded mapreduce

GeometryOpsCore/src/types/applicators.jl

@@ -54,8 +54,26 @@ const WithXYZM = ApplyToPoint{true,true}
 
 (t::WithXY)(p) = t.f(GI.x(p), GI.y(p))
 (t::WithXYZ)(p) = t.f(GI.x(p), GI.y(p), GI.z(p))
-(t::WithXYZM)(p) = t.f(GI.x(p), GI.y(p), GI.m(p))
-(t::WithXYM)(p) = t.f(GI.x(p), GI.y(p), GI.z(p), GI.m(p))
+(t::WithXYZM)(p) = t.f(GI.x(p), GI.y(p), GI.z(p), GI.m(p))
+(t::WithXYM)(p) = t.f(GI.x(p), GI.y(p), GI.m(p))
+
+"""
+    WithTrait(f)
+
+WithTrait is a functor that applies a function to a trait and an object.
+
+Specifically, the calling convention is for `f` is changed
+from `f(geom)` to `f(trait, geom; kw...)`.
+
+This is useful to keep the trait materialized through the call stack,
+which can improve inferrability and performance.
+"""
+struct WithTrait{F} <: Applicator{F, Nothing}
+    f::F
+end
+
+(a::WithTrait)(trait::GI.AbstractTrait, obj; kw...) = a.f(trait, obj; kw...)
+rebuild(::WithTrait, f::F) where {F} = WithTrait{F}(f)
 
 # ***
 

Project.toml

@@ -1,7 +1,7 @@
 name = "GeometryOps"
 uuid = "3251bfac-6a57-4b6d-aa61-ac1fef2975ab"
 authors = ["Anshul Singhvi <[email protected]>", "Rafael Schouten <[email protected]>", "Skylar Gering <[email protected]>", "and contributors"]
-version = "0.1.23"
+version = "0.1.24"
 
 [deps]
 AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -13,20 +13,22 @@ ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110"
 Extents = "411431e0-e8b7-467b-b5e0-f676ba4f2910"
 GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
 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"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
 [weakdeps]
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 FlexiJoins = "e37f2e79-19fa-4eb7-8510-b63b51fe0a37"
 LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 TGGeometry = "d7e755d2-3c95-4bcf-9b3c-79ab1a78647b"
 
 [extensions]
+GeometryOpsDataFramesExt = "DataFrames"
 GeometryOpsFlexiJoinsExt = "FlexiJoins"
 GeometryOpsLibGEOSExt = "LibGEOS"
 GeometryOpsProjExt = "Proj"
@@ -37,14 +39,15 @@ AbstractTrees = "0.4"
 AdaptivePredicates = "1.2"
 CoordinateTransformations = "0.5, 0.6"
 DataAPI = "1"
+DataFrames = "1"
 DelaunayTriangulation = "1.0.4"
 ExactPredicates = "2.2.8"
 Extents = "0.1.5"
 FlexiJoins = "0.1.30"
 GeoFormatTypes = "0.4"
 GeoInterface = "1.2"
 GeometryBasics = "0.4.7, 0.5"
-GeometryOpsCore = "=0.1.5"
+GeometryOpsCore = "=0.1.6"
 LibGEOS = "0.9.2"
 LinearAlgebra = "1"
 Proj = "1"
@@ -63,6 +66,7 @@ Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 FlexiJoins = "e37f2e79-19fa-4eb7-8510-b63b51fe0a37"
 GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9"
+GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
 NaturalEarth = "436b0209-26ab-4e65-94a9-6526d86fea76"
@@ -77,4 +81,4 @@ TGGeometry = "d7e755d2-3c95-4bcf-9b3c-79ab1a78647b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["ArchGDAL", "CoordinateTransformations", "DataFrames", "Distributions", "DimensionalData", "Downloads", "FlexiJoins", "GeoJSON", "Proj", "JLD2", "LibGEOS", "Random", "Rasters", "NaturalEarth", "OffsetArrays", "Polylabel", "SafeTestsets", "Shapefile", "TGGeometry", "Test"]
+test = ["ArchGDAL", "CoordinateTransformations", "DataFrames", "Distributions", "DimensionalData", "Downloads", "FlexiJoins", "GeoJSON", "GeometryBasics", "Proj", "JLD2", "LibGEOS", "Random", "Rasters", "NaturalEarth", "OffsetArrays", "Polylabel", "SafeTestsets", "Shapefile", "TGGeometry", "Test"]

benchmarks/benchmarks.jl

@@ -62,87 +62,6 @@ end
 plot_trials(circle_area_suite)
 
 
-# ## Vancouver watershed benchmarks
-#=
-
-Vancouver Island has ~1,300 watersheds.  LibGEOS uses this exact data
-in their tests, so we'll use it in ours as well!
-
-We'll start by loading the data, and then we'll use it to benchmark various operations.
-
-=#
-
-# The CRS for this file is EPSG:3005, or as a PROJ string,
-# `"+proj=aea +lat_1=50 +lat_2=58.5 +lat_0=45 +lon_0=-126 +x_0=1000000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"`
-# TODO: this doesn't work with WellKnownGeometry.  Why?
-wkt_geoms = LG.readgeom.(eachline("/Users/anshul/Downloads/watersheds.wkt"), (LG.WKTReader(LG.get_global_context()),))
-vancouver_polygons = GI.getgeom.(wkt_geoms, 1); #.|> GO.tuples;
-
-import SortTileRecursiveTree as STR
-tree = STR.STRtree(vancouver_polygons)
-query_result = STR.query(tree, GI.extent(vancouver_polygons[1]))
-
-GO.intersects.((vancouver_polygons[1],), vancouver_polygons[query_result])
-
-go_vp = GO.tuples.(vancouver_polygons[1:2])
-@be GO.union($(go_vp[1]), $(go_vp[2]); target = $GI.PolygonTrait())
-@be LG.union($(vancouver_polygons[1]), $(vancouver_polygons[2]))
-
-all_intersected = falses(length(vancouver_polygons))
-accumulator = deepcopy(vancouver_polygons[1])
-all_intersected[1] = true
-i = 1
-# query_result = STR.query(tree, GI.extent(accumulator))
-# for idx in query_result
-#     println("Assessing $idx")
-#     if !all_intersected[idx] && LG.intersects(vancouver_polygons[idx], accumulator)
-#         println("Assimilating $idx")
-#         result = LG.union(vancouver_polygons[idx], accumulator#=; target = GI.PolygonTrait()=#)
-#         # @show length(result)
-#         accumulator = result#[1]
-#         all_intersected[idx] = true
-#     end
-# end
-display(poly(vancouver_polygons[all_intersected]; color = rand(RGBf, sum(all_intersected))))
-display(poly(accumulator))
-@time while !(all(all_intersected)) && i < length(vancouver_polygons)
-    println("Began iteration $i")
-    query_result = STR.query(tree, GI.extent(accumulator))
-    for idx in query_result
-        if !(all_intersected[idx] || !(LG.intersects(vancouver_polygons[idx], accumulator)))
-            println("Assimilating $idx")
-            result = LG.union(vancouver_polygons[idx], accumulator#=; target = GI.PolygonTrait()=#)
-            # @show length(result)
-            accumulator = result#[1]
-            all_intersected[idx] = true
-        end
-    end
-    display(poly(vancouver_polygons[all_intersected]; color = rand(RGBf, sum(all_intersected))))
-    println("Finished iteration $i")
-    # wait_for_key("Press any key to continue to the next iteration.")
-    i += 1
-end 
-
-fig = Figure()
-ax = Axis(fig[1, 1]; title = "STRTree query for polygons", aspect = DataAspect())
-for (i, result_index) in enumerate(query_result)
-    poly!(ax, vancouver_polygons[result_index]; color = Makie.wong_colors()[i], label = "$result_index")
-end
-Legend(fig[1, 2], ax)
-fig
-
-
-
-# TODO: 
-# - Vancouver watersheds:
-#    - Intersection on pre-buffered geometry
-#    - Polygon union by reduction (perhaps pre-sort by border order, so we don't end up with useless polygons)
-#    - Queries using STRTree.jl
-#    - Potentially using a prepared geometry based approach to multithreaded reductive joining
-#    - Implement multipolygon joining.  How?  Query intersection or touching for each individual geometry,
-#      and implement a 
-
-
 ## Segmentization
 
 
@@ -180,9 +99,9 @@ n_points_values = round.(Int, exp10.(LinRange(0.7, 6, 15)))
     circle_union_suite["LibGEOS"][n_points]     = @be LG.union($lg_circle_left, $lg_circle_right)
 end
 
-plot_trials(circle_difference_suite)
-plot_trials(circle_intersection_suite)
-plot_trials(circle_union_suite)
+plot_trials(circle_difference_suite; legend_position = (2, 1), legend_kws = (; orientation = :horizontal, nbanks = 2))
+plot_trials(circle_intersection_suite; legend_position = (2, 1), legend_kws = (; orientation = :horizontal, nbanks = 2))
+plot_trials(circle_union_suite; legend_position = (2, 1), legend_kws = (; orientation = :horizontal, nbanks = 2))
 
 usa_poly_suite = BenchmarkGroup()
 usa_difference_suite = usa_poly_suite["difference"] = BenchmarkGroup(["title:USA difference", "subtitle:Tested on CONUS"])

docs/Project.toml

@@ -36,6 +36,7 @@ Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 MakieThemes = "e296ed71-da82-5faf-88ab-0034a9761098"
 Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
 MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca"
+Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 MultiFloats = "bdf0d083-296b-4888-a5b6-7498122e68a5"
 NaturalEarth = "436b0209-26ab-4e65-94a9-6526d86fea76"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"