Add support for tables to apply and applyreduce (#91)

* Add support for tables to `apply`

* Fix the `reproject` implementation to use constructed traits

* Fix error when testing `threaded`

* Comment, fix indentation in apply

* Implement applyreduce

* document in readme and docs that tables are supported

* Add DataFrames and Shapefile to project and docs/project

* Remove unused function

Co-authored-by: Rafael Schouten <[email protected]>

* Add some truly wacko tests

* Fix test

* `_to_table` -> `_table`

* Annotate that `_is3d` only returns bools

Apparently this was not obvious to the compiler

* Don't be lazy when getting names

this promotes type stability and I guess we don't care about allocs anyway.

---------

Co-authored-by: Rafael Schouten <[email protected]>

Author: asinghvi17

Project.toml

@@ -9,6 +9,7 @@ GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
 [weakdeps]
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
@@ -23,19 +24,22 @@ GeometryBasics = "0.4.7"
 LinearAlgebra = "1"
 Proj = "1"
 Statistics = "1"
+Tables = "1"
 julia = "1.9"
 
 [extras]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
 CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
 GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
 Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["ArchGDAL", "CoordinateTransformations", "Distributions", "GeoFormatTypes", "GeoJSON", "Proj", "JLD2", "LibGEOS", "Random", "Test"]
+test = ["ArchGDAL", "CoordinateTransformations", "DataFrames", "Distributions", "GeoFormatTypes", "GeoJSON", "Proj", "JLD2", "LibGEOS", "Random", "Shapefile", "Test"]

README.md

@@ -24,18 +24,19 @@ GeometryOps tries to offer most of the basic geometry operations you'd need, imp
 
 - General geometry methods (OGC methods): `equals`, `extent`, `distance`, `crosses`, `contains`, `intersects`, etc
 - Targeted function application over large nested geometries (`apply`) and reduction over geometries (`applyreduce`)
+    - Both `apply` and `applyreduce` consume arbitrary tables as well, like DataFrames!
 - `signed_area`, `centroid`, `distance`, etc for valid geometries
 - Line and polygon simplification (`simplify`)
 - Polygon clipping, `intersection`, `difference` and `union`
 - Generalized barycentric coordinates in polygons (`barycentric_coordinates`)
 - Projection of geometries between coordinate reference systems using [Proj.jl](https://github.com/JuliaGeo/Proj.jl)
 - Polygonization of raster images by contour detection (`polygonize`)
+- Segmentization/densification of geometry, both linearly and by geodesic paths (`segmentize`)
 
 See the "API" page in the docs for a more complete list!
 
 ### Planned additions
 
-- Arclength interpolation (absolute and relative)
 - Buffering, hulls (convex and otherwise)
 - Checks for valid geometries (empty linestrings, null points, etc) ([#14](https://github.com/asinghvi17/GeometryOps.jl/issues/14))
 - Operations on spherical (non-Euclidean) geometry ([#17](https://github.com/asinghvi17/GeometryOps.jl/issues/17))

docs/Project.toml

@@ -4,6 +4,7 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Chairmarks = "0ca39b1e-fe0b-4e98-acfc-b1656634c4de"
 CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DataInterpolations = "82cc6244-b520-54b8-b5a6-8a565e85f1d0"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
@@ -17,8 +18,9 @@ GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
 GeometryOps = "3251bfac-6a57-4b6d-aa61-ac1fef2975ab"
 LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
-Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 MakieThemes = "e296ed71-da82-5faf-88ab-0034a9761098"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
+Proj = "c94c279d-25a6-4763-9509-64d165bea63e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"

docs/src/introduction.md

@@ -20,6 +20,8 @@ The `apply` function allows you to decompose a given collection of geometries do
 
 Functionally, it's similar to `map` in the way you apply it to geometries.
 
+`apply` and `applyreduce` take any geometry, vector of geometries, collection of geometries, or table (like `Shapefile.Table`, `DataFrame`, or `GeoTable`)!
+
 ### What's this `GeoInterface.Wrapper` thing?
 
 Write a comment about GeoInterface.Wrapper and why it helps in type stability to guarantee a particular return type.

ext/GeometryOpsProjExt/reproject.jl

@@ -39,11 +39,11 @@ function reproject(geom, source_crs, target_crs;
 end
 function reproject(geom, transform::Proj.Transformation; time=Inf, target_crs=nothing, kw...)
     if _is3d(geom)
-        return apply(GI.PointTrait, geom; crs=target_crs, kw...) do p
+        return apply(GI.PointTrait(), geom; crs=target_crs, kw...) do p
             transform(GI.x(p), GI.y(p), GI.z(p))
         end
     else
-        return apply(GI.PointTrait, geom; crs=target_crs, kw...) do p
+        return apply(GI.PointTrait(), geom; crs=target_crs, kw...) do p
             transform(GI.x(p), GI.y(p))
         end
     end

src/GeometryOps.jl

@@ -4,6 +4,7 @@ module GeometryOps
 
 using GeoInterface
 using GeometryBasics
+import Tables
 using LinearAlgebra, Statistics
 import GeometryBasics.StaticArrays
 import Base.@kwdef

src/primitives.jl

@@ -212,17 +212,51 @@ 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; threaded, kw...) where F
-    if threaded
-        # `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)
+@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.
+    return Tables.materializer(iterable)(
+        merge(
+            NamedTuple{(geometry_column,), Base.Tuple{typeof(new_geometry)}}((new_geometry,)),
+            NamedTuple(Iterators.map(_get_col_pair, new_names))
+        )
+    )
+end
+
 # Rewrap all FeatureCollectionTrait feature collections as GI.FeatureCollection
-# Maybe use threads to call _apply on componenet features
+# 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
@@ -325,14 +359,34 @@ end
     _mapreducetasks(applyreduce_array, op, eachindex(A), threaded; init)
 end
 # Try to applyreduce over iterables
-@inline function _applyreduce(f::F, op::O, target, ::Nothing, iterable; threaded, init) where {F, O}
-    applyreduce_iterable(i) = _applyreduce(f, op, target, x; threaded=_False(), init)
-    if threaded # Try to `collect` and reduce over the vector with threads
-        _applyreduce(f, op, target, collect(iterable); threaded, init)
+@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
-        # Try to `mapreduce` the iterable as-is
-        mapreduce(applyreduce_iterable, op, iterable; init)
+        applyreduce_iterable(i) = _applyreduce(f, op, target, x; 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}

src/utils.jl

@@ -1,10 +1,10 @@
 # # Utility functions
 
-_is3d(geom) = _is3d(GI.trait(geom), geom)
-_is3d(::GI.AbstractGeometryTrait, geom) = GI.is3d(geom)
-_is3d(::GI.FeatureTrait, feature) = _is3d(GI.geometry(feature))
-_is3d(::GI.FeatureCollectionTrait, fc) = _is3d(GI.getfeature(fc, 1))
-_is3d(::Nothing, geom) = _is3d(first(geom)) # Otherwise step into an itererable
+_is3d(geom)::Bool = _is3d(GI.trait(geom), geom)
+_is3d(::GI.AbstractGeometryTrait, geom)::Bool = GI.is3d(geom)
+_is3d(::GI.FeatureTrait, feature)::Bool = _is3d(GI.geometry(feature))
+_is3d(::GI.FeatureCollectionTrait, fc)::Bool = _is3d(GI.getfeature(fc, 1))
+_is3d(::Nothing, geom)::Bool = _is3d(first(geom)) # Otherwise step into an itererable
 
 _npoint(x) = _npoint(trait(x), x)
 _npoint(::Nothing, xs::AbstractArray) = sum(_npoint, xs)

test/primitives.jl

@@ -4,6 +4,7 @@ import GeoInterface as GI
 import GeometryOps as GO
 import GeometryBasics as GB
 import Proj
+import Shapefile, DataFrames
 
 pv1 = [(1, 2), (3, 4), (5, 6), (1, 2)]
 pv2 = [(3, 4), (5, 6), (6, 7), (3, 4)]
@@ -19,8 +20,49 @@ poly = GI.Polygon([lr1, lr2])
 
     @test flipped_poly == GI.Polygon([GI.LinearRing([(2, 1), (4, 3), (6, 5), (2, 1)]), 
                                       GI.LinearRing([(4, 3), (6, 5), (7, 6), (4, 3)])])
+
+    @testset "Tables.jl support" begin
+        mktempdir() do dir
+        cd(dir) do
+
+            download("https://rawcdn.githack.com/nvkelso/natural-earth-vector/v5.1.2/110m_cultural/ne_110m_admin_0_countries.shp", "countries.shp")
+            download("https://rawcdn.githack.com/nvkelso/natural-earth-vector/v5.1.2/110m_cultural/ne_110m_admin_0_countries.shx", "countries.shx")
+            download("https://rawcdn.githack.com/nvkelso/natural-earth-vector/v5.1.2/110m_cultural/ne_110m_admin_0_countries.dbf", "countries.dbf")
+            download("https://rawcdn.githack.com/nvkelso/natural-earth-vector/v5.1.2/110m_cultural/ne_110m_admin_0_countries.prj", "countries.prj")
+            countries_table = Shapefile.Table("countries.shp")
+
+            @testset "Shapefile" begin
+                centroid_table = GO.apply(GO.centroid, GO.TraitTarget(GI.PolygonTrait(), GI.MultiPolygonTrait()), countries_table);
+                centroid_geometry = centroid_table.geometry
+                # Test that the centroids are correct
+                @test all(centroid_geometry .== GO.centroid.(countries_table.geometry))
+                @testset "Columns are preserved" begin  
+                    for column in Iterators.filter(!=(:geometry), GO.Tables.columnnames(countries_table))
+                        @test all(GO.Tables.getcolumn(centroid_table, column) .== GO.Tables.getcolumn(countries_table, column))
+                    end
+                end
+            end
+
+            @testset "DataFrames" begin
+                countries_df = DataFrames.DataFrame(countries_table)
+                centroid_df = GO.apply(GO.centroid, GO.TraitTarget(GI.PolygonTrait(), GI.MultiPolygonTrait()), countries_df);
+                @test centroid_df isa DataFrames.DataFrame
+                centroid_geometry = centroid_df.geometry
+                # Test that the centroids are correct
+                @test all(centroid_geometry .== GO.centroid.(countries_df.geometry))
+                @testset "Columns are preserved" begin  
+                    for column in Iterators.filter(!=(:geometry), GO.Tables.columnnames(countries_df))
+                        @test all(centroid_df[!, column] .== countries_df[!, column])
+                    end
+                end
+            end
+        end
+        end
+    end
 end
 
+
+
 @testset "unwrap" begin
     flipped_vectors = GO.unwrap(GI.PointTrait, poly) do p
         (GI.y(p), GI.x(p))