Improve apply on tables, add a one-crs method to reproject

GeometryOpsCore/src/apply.jl

@@ -175,20 +175,40 @@ 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}
+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.
-    geometry_column = first(GI.geometrycolumns(iterable))
-    new_geometry = _apply(f, target, Tables.getcolumn(iterable, geometry_column); threaded, kw...)
+    geometry_columns = if isnothing(geometrycolumn)
+        GI.geometrycolumns(iterable)
+    elseif geometrycolumn isa NTuple{N, <: Symbol} where N
+        geometrycolumn
+    elseif geometrycolumn isa Symbol
+        (geometrycolumn,)
+    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)
+        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))
+            """
+            ))
+    end
+    new_geometry_vecs = map(geometry_columns) do colname
+        _apply(f, target, Tables.getcolumn(iterable, colname); threaded, kw...)
+    end
     # 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)
+    new_names = filter(x -> !(x in geometry_columns), 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{geometry_columns, Base.Tuple{typeof.(new_geometry_vecs)...}}(new_geometry_vecs),
             NamedTuple(Iterators.map(_get_col_pair, new_names))
         )
     )
@@ -201,7 +221,9 @@ function _apply_table(f::F, target, iterable::IterableType; threaded, kw...) whe
         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
+                if !preserve_default_metadata 
+                    style == :default && continue
+                end
                 # We assume that any other style is preserved.
                 DataAPI.metadata!(result, key, value; style)
             end
@@ -214,16 +236,11 @@ function _apply_table(f::F, target, iterable::IterableType; threaded, kw...) whe
         # 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)
+            DataAPI.metadata!(result, "GEOINTERFACE:geometrycolumns", geometry_columns; style = :note)
         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)
+        new_crs = get(kw, :crs, GI.crs(iterable)) 
+        DataAPI.metadata!(result, "GEOINTERFACE:crs", new_crs; style = :note)
     end
 
     return result

ext/GeometryOpsProjExt/reproject.jl

@@ -12,15 +12,20 @@
 )
     if isnothing(transform)
         if isnothing(source_crs) 
-            source_crs = if GI.trait(geom) isa Nothing && geom isa AbstractArray
-                GeoInterface.crs(first(geom))
-            else
-                GeoInterface.crs(geom)
+            # this will check DataAPI.jl metadata as well
+            source_crs = GI.crs(geom)
+            # if GeoInterface somehow missed the CRS, we assume it can only
+            # be an iterable, because GeoInterface queries DataAPI.jl metadata
+            # from tables and such things.
+            if isnothing(source_crs) && isnothing(GI.trait(geom))
+                if Base.isiterable(typeof(geom))
+                    source_crs = GI.crs(first(geom))
+                end
             end
         end
 
         # If its still nothing, error
-        isnothing(source_crs) && throw(ArgumentError("geom has no crs attached. Pass a `source_crs` keyword"))
+        isnothing(source_crs) && throw(ArgumentError("geom has no crs attached. Pass a `source_crs` keyword argument."))
 
         # Otherwise reproject
         reproject(geom, source_crs, target_crs; kw...)
@@ -34,10 +39,27 @@
 end
 function reproject(
     geom, source_crs::CRSType, target_crs::CRSType; always_xy=true, kw...
-) where CRSType <: Union{GeoFormatTypes.GeoFormat, Proj.CRS}
+) where CRSType <: Union{GeoFormatTypes.GeoFormat, Proj.CRS, String}
     transform = Proj.Transformation(source_crs, target_crs; always_xy)
     return reproject(geom, transform; target_crs, kw...)
 end
+function reproject(
+    geom, target_crs::CRSType; kw...
+) where CRSType <: Union{GeoFormatTypes.GeoFormat, Proj.CRS, String}
+    source_crs = GI.crs(geom)
+    if isnothing(source_crs) 
+        if GI.DataAPI.metadatasupport(typeof(geom)).read
+            source_crs = GI.crs(geom)
+        end
+        if isnothing(source_crs)
+            if geom isa AbstractArray
+                source_crs = GI.crs(first(geom))
+            end
+        end
+    end
+    isnothing(source_crs) && throw(ArgumentError("geom has no crs attached. Pass a `source_crs` before the current target crs you have passed."))
+    return reproject(geom; source_crs, target_crs, kw...)
+end
 function reproject(geom, transform::Proj.Transformation; 
     context=C_NULL, 
     target_crs=nothing, 

src/utils/utils.jl

@@ -1,10 +1,31 @@
 # # Utility functions
 
-_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
+_is3d(geom; kw...)::Bool = _is3d(GI.trait(geom), geom; kw...)
+_is3d(::GI.AbstractGeometryTrait, geom; geometrycolumn = nothing)::Bool = GI.is3d(geom)
+_is3d(::GI.FeatureTrait, feature; geometrycolumn = nothing)::Bool = _is3d(GI.geometry(feature))
+_is3d(::GI.FeatureCollectionTrait, fc; geometrycolumn = nothing)::Bool = _is3d(GI.getfeature(fc, 1))
+function _is3d(::Nothing, geom; geometrycolumn = nothing)::Bool
+    if Tables.istable(geom)
+        geometrycolumn = isnothing(geometrycolumn) ? GI.geometrycolumns(geom) : geometrycolumn isa Symbol ? (geometrycolumn,) : geometrycolumn
+        # take the first geometry column
+        # TODO: this is a bad guess - this should really be on the vector level somehow.  
+        # Maybe a configurable applicator again....
+        first_geom = if Tables.rowaccess(geom)
+            Tables.getcolumn(first(Tables.rows(geom)), first(geometrycolumn))
+        else # column access assumed
+            first(Tables.getcolumn(geom, first(geometrycolumn)))
+        end
+        return _is3d(first_geom)
+    else # assume iterable
+        first_geom = first(geom)
+        if GI.trait(first_geom) isa GI.AbstractTrait
+            return _is3d(first_geom)
+        else
+            return false # couldn't figure it out!
+        end
+    end
+
+end
 
 _npoint(x) = _npoint(trait(x), x)
 _npoint(::Nothing, xs::AbstractArray) = sum(_npoint, xs)

test/primitives.jl

@@ -56,30 +56,45 @@ poly = GI.Polygon([lr1, lr2])
                     end
                 end
 
-                @testset "DataFrames" begin
-                    countries_df = DataFrames.DataFrame(countries_table)
-                    GO.DataAPI.metadata!(countries_df, "note metadata", "note metadata value"; style = :note)
-                    GO.DataAPI.metadata!(countries_df, "default metadata", "default metadata value"; style = :default)
-                    centroid_df = GO.apply(GO.centroid, GO.TraitTarget(GI.PolygonTrait(), GI.MultiPolygonTrait()), countries_df; crs = GFT.EPSG(3031));
-                    @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(missing_or_equal.(centroid_df[!, column], countries_df[!, column]))
-                        end
-                    end
-                    @testset "Metadata preservation (or not)" begin
-                        @test DataAPI.metadata(centroid_df, "note metadata") == "note metadata value"
-                        @test !("default metadata" in DataAPI.metadatakeys(centroid_df))
-                        @test DataAPI.metadata(centroid_df, "GEOINTERFACE:geometrycolumns") == (:geometry,)
-                        @test DataAPI.metadata(centroid_df, "GEOINTERFACE:crs") == GFT.EPSG(3031)
+            @testset "DataFrames" begin
+                countries_df = DataFrames.DataFrame(countries_table)
+                GO.DataAPI.metadata!(countries_df, "note metadata", "note metadata value"; style = :note)
+                GO.DataAPI.metadata!(countries_df, "default metadata", "default metadata value"; style = :default)
+                centroid_df = GO.apply(GO.centroid, GO.TraitTarget(GI.PolygonTrait(), GI.MultiPolygonTrait()), countries_df; crs = GFT.EPSG(3031));
+                # Test that the Tables.jl materializer is used
+                @test centroid_df isa DataFrames.DataFrame
+                # Test that the centroids are correct
+                @test all(centroid_df.geometry .== GO.centroid.(countries_df.geometry))
+                @testset "Columns are preserved" begin  
+                    for column in filter(!=(:geometry), GO.Tables.columnnames(countries_df))
+                        @test all(missing_or_equal.(centroid_df[!, column], countries_df[!, column]))
                     end
                 end
+                @testset "Multiple geometry columns in metadata" begin
+                    # set up a dataframe with multiple geometry columns
+                    countries_df2 = deepcopy(countries_df)
+                    countries_df2.centroid = GO.centroid.(countries_df2.geometry)
+                    GI.DataAPI.metadata!(countries_df2, "GEOINTERFACE:geometrycolumns", (:geometry, :centroid); style = :note)
+                    transformed = GO.transform(p -> p .+ 3, countries_df2)
+                    @test GI.DataAPI.metadata(transformed, "GEOINTERFACE:geometrycolumns") == (:geometry, :centroid)
+                    @test GI.DataAPI.metadata(transformed, "GEOINTERFACE:crs") == GI.crs(countries_df2)
+                    # Test that the transformation was actually applied to both geometry columns.
+                    @test all(map(zip(countries_df2.geometry, transformed.geometry)) do (o, n)
+                        GO.equals(GO.transform(p -> p .+ 3, o), n)
+                    end)
+                    @test all(map(zip(countries_df2.centroid, transformed.centroid)) do (o, n)
+                        any(isnan, o) || GO.equals(GO.transform(p -> p .+ 3, o), n)
+                    end)
+                end
+                end
             end
         end
         end
+        @testset "Wrong geometry column kwarg" begin
+            tab = Tables.dictcolumntable((; geometry = [(1, 2), (3, 4), (5, 6)], other = [1, 2, 3]))
+            @test_throws "got a Float64" GO.transform(identity, tab; geometrycolumn = 1000.0)
+            @test_throws "but the table has columns" GO.transform(identity, tab; geometrycolumn = :somethingelse)
+        end
     end
 end
 

test/transformations/reproject.jl

@@ -79,3 +79,16 @@ _xy(p) = GI.x(p), GI.y(p)
 
     GO.reproject(multipolygon4326; target_crs=ProjString("+proj=moll +type=crs"))
 end
+
+
+@testset "Reproject with only target crs" begin
+    multipolygon4326 = GO.tuples(multipolygon; crs = EPSG(4326))
+    multipolygon3857 = GO.reproject(multipolygon, EPSG(4326), EPSG(3857))
+    @test GI.crs(GO.reproject(multipolygon3857; target_crs=EPSG(4326))) == EPSG(4326)
+    @test GI.crs(GO.reproject(multipolygon3857, EPSG(4326))) == EPSG(4326)
+
+    v = [multipolygon3857]
+    @test GI.crs(only(GO.reproject(v, EPSG(4326)))) == EPSG(4326)
+end
+
+