Skip to content

Commit 7f82f02

Browse files
asinghvi17rafaqz
asinghvi17
and
rafaqz
authored
Improve apply on tables, add a one-crs method to reproject (#302)
Co-authored-by: Rafael Schouten <[email protected]>
1 parent e7e5240 commit 7f82f02

File tree

5 files changed

+132
-44
lines changed

5 files changed

+132
-44
lines changed

‎GeometryOpsCore/src/apply.jl

Lines changed: 31 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -175,20 +175,40 @@ with the same schema, but with the new geometry column.
175175
This new table may be of the same type as the old one iff `Tables.materializer` is defined for
176176
that table. If not, then a `NamedTuple` is returned.
177177
=#
178-
function _apply_table(f::F, target, iterable::IterableType; threaded, kw...) where {F, IterableType}
178+
function _apply_table(f::F, target, iterable::IterableType; geometrycolumn = nothing, preserve_default_metadata = false, threaded, kw...) where {F, IterableType}
179179
_get_col_pair(colname) = colname => Tables.getcolumn(iterable, colname)
180180
# We extract the geometry column and run `apply` on it.
181-
geometry_column = first(GI.geometrycolumns(iterable))
182-
new_geometry = _apply(f, target, Tables.getcolumn(iterable, geometry_column); threaded, kw...)
181+
geometry_columns = if isnothing(geometrycolumn)
182+
GI.geometrycolumns(iterable)
183+
elseif geometrycolumn isa NTuple{N, <: Symbol} where N
184+
geometrycolumn
185+
elseif geometrycolumn isa Symbol
186+
(geometrycolumn,)
187+
else
188+
throw(ArgumentError("geometrycolumn must be a Symbol or a tuple of Symbols, got a $(typeof(geometrycolumn))"))
189+
end
190+
if !all(Base.Fix2(in, Tables.columnnames(iterable)), geometry_columns)
191+
throw(ArgumentError(
192+
"""
193+
`apply`: the `geometrycolumn` kwarg must be a subset of the column names of the table,
194+
got $(geometry_columns)
195+
but the table has columns
196+
$(Tables.columnnames(iterable))
197+
"""
198+
))
199+
end
200+
new_geometry_vecs = map(geometry_columns) do colname
201+
_apply(f, target, Tables.getcolumn(iterable, colname); threaded, kw...)
202+
end
183203
# Then, we obtain the schema of the table,
184204
old_schema = Tables.schema(iterable)
185205
# filter the geometry column out,
186-
new_names = filter(Base.Fix1(!==, geometry_column), old_schema.names)
206+
new_names = filter(x -> !(x in geometry_columns), old_schema.names)
187207
# and try to rebuild the same table as the best type - either the original type of `iterable`,
188208
# or a named tuple which is the default fallback.
189209
result = Tables.materializer(iterable)(
190210
merge(
191-
NamedTuple{(geometry_column,), Base.Tuple{typeof(new_geometry)}}((new_geometry,)),
211+
NamedTuple{geometry_columns, Base.Tuple{typeof.(new_geometry_vecs)...}}(new_geometry_vecs),
192212
NamedTuple(Iterators.map(_get_col_pair, new_names))
193213
)
194214
)
@@ -201,7 +221,9 @@ function _apply_table(f::F, target, iterable::IterableType; threaded, kw...) whe
201221
if DataAPI.metadatasupport(IterableType).read
202222
for (key, (value, style)) in DataAPI.metadata(iterable; style = true)
203223
# Default styles are not preserved on data transformation, so we must skip them!
204-
style == :default && continue
224+
if !preserve_default_metadata
225+
style == :default && continue
226+
end
205227
# We assume that any other style is preserved.
206228
DataAPI.metadata!(result, key, value; style)
207229
end
@@ -214,16 +236,11 @@ function _apply_table(f::F, target, iterable::IterableType; threaded, kw...) whe
214236
# so we don't need to set them.
215237
if !("GEOINTERFACE:geometrycolumns" in mdk)
216238
# If the geometry columns are not already set, we need to set them.
217-
DataAPI.metadata!(result, "GEOINTERFACE:geometrycolumns", (geometry_column,); style = :default)
239+
DataAPI.metadata!(result, "GEOINTERFACE:geometrycolumns", geometry_columns; style = :note)
218240
end
219241
# Force reset CRS always, since you can pass `crs` to `apply`.
220-
new_crs = if haskey(kw, :crs)
221-
kw[:crs]
222-
else
223-
GI.crs(iterable) # this will automatically check `GEOINTERFACE:crs` unless the type has a specialized implementation.
224-
end
225-
226-
DataAPI.metadata!(result, "GEOINTERFACE:crs", new_crs; style = :default)
242+
new_crs = get(kw, :crs, GI.crs(iterable))
243+
DataAPI.metadata!(result, "GEOINTERFACE:crs", new_crs; style = :note)
227244
end
228245

229246
return result

‎ext/GeometryOpsProjExt/reproject.jl

Lines changed: 28 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -12,15 +12,20 @@ function reproject(geom;
1212
)
1313
if isnothing(transform)
1414
if isnothing(source_crs)
15-
source_crs = if GI.trait(geom) isa Nothing && geom isa AbstractArray
16-
GeoInterface.crs(first(geom))
17-
else
18-
GeoInterface.crs(geom)
15+
# this will check DataAPI.jl metadata as well
16+
source_crs = GI.crs(geom)
17+
# if GeoInterface somehow missed the CRS, we assume it can only
18+
# be an iterable, because GeoInterface queries DataAPI.jl metadata
19+
# from tables and such things.
20+
if isnothing(source_crs) && isnothing(GI.trait(geom))
21+
if Base.isiterable(typeof(geom))
22+
source_crs = GI.crs(first(geom))
23+
end
1924
end
2025
end
2126

2227
# If its still nothing, error
23-
isnothing(source_crs) && throw(ArgumentError("geom has no crs attached. Pass a `source_crs` keyword"))
28+
isnothing(source_crs) && throw(ArgumentError("geom has no crs attached. Pass a `source_crs` keyword argument."))
2429

2530
# Otherwise reproject
2631
reproject(geom, source_crs, target_crs; kw...)
@@ -34,10 +39,27 @@ function reproject(geom, source_crs, target_crs; always_xy=true, kw...)
3439
end
3540
function reproject(
3641
geom, source_crs::CRSType, target_crs::CRSType; always_xy=true, kw...
37-
) where CRSType <: Union{GeoFormatTypes.GeoFormat, Proj.CRS}
42+
) where CRSType <: Union{GeoFormatTypes.GeoFormat, Proj.CRS, String}
3843
transform = Proj.Transformation(source_crs, target_crs; always_xy)
3944
return reproject(geom, transform; target_crs, kw...)
4045
end
46+
function reproject(
47+
geom, target_crs::CRSType; kw...
48+
) where CRSType <: Union{GeoFormatTypes.GeoFormat, Proj.CRS, String}
49+
source_crs = GI.crs(geom)
50+
if isnothing(source_crs)
51+
if GI.DataAPI.metadatasupport(typeof(geom)).read
52+
source_crs = GI.crs(geom)
53+
end
54+
if isnothing(source_crs)
55+
if geom isa AbstractArray
56+
source_crs = GI.crs(first(geom))
57+
end
58+
end
59+
end
60+
isnothing(source_crs) && throw(ArgumentError("geom has no crs attached. Pass a `source_crs` before the current target crs you have passed."))
61+
return reproject(geom; source_crs, target_crs, kw...)
62+
end
4163
function reproject(geom, transform::Proj.Transformation;
4264
context=C_NULL,
4365
target_crs=nothing,

‎src/utils/utils.jl

Lines changed: 26 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,10 +1,31 @@
11
# # Utility functions
22

3-
_is3d(geom)::Bool = _is3d(GI.trait(geom), geom)
4-
_is3d(::GI.AbstractGeometryTrait, geom)::Bool = GI.is3d(geom)
5-
_is3d(::GI.FeatureTrait, feature)::Bool = _is3d(GI.geometry(feature))
6-
_is3d(::GI.FeatureCollectionTrait, fc)::Bool = _is3d(GI.getfeature(fc, 1))
7-
_is3d(::Nothing, geom)::Bool = _is3d(first(geom)) # Otherwise step into an itererable
3+
_is3d(geom; kw...)::Bool = _is3d(GI.trait(geom), geom; kw...)
4+
_is3d(::GI.AbstractGeometryTrait, geom; geometrycolumn = nothing)::Bool = GI.is3d(geom)
5+
_is3d(::GI.FeatureTrait, feature; geometrycolumn = nothing)::Bool = _is3d(GI.geometry(feature))
6+
_is3d(::GI.FeatureCollectionTrait, fc; geometrycolumn = nothing)::Bool = _is3d(GI.getfeature(fc, 1))
7+
function _is3d(::Nothing, geom; geometrycolumn = nothing)::Bool
8+
if Tables.istable(geom)
9+
geometrycolumn = isnothing(geometrycolumn) ? GI.geometrycolumns(geom) : geometrycolumn isa Symbol ? (geometrycolumn,) : geometrycolumn
10+
# take the first geometry column
11+
# TODO: this is a bad guess - this should really be on the vector level somehow.
12+
# Maybe a configurable applicator again....
13+
first_geom = if Tables.rowaccess(geom)
14+
Tables.getcolumn(first(Tables.rows(geom)), first(geometrycolumn))
15+
else # column access assumed
16+
first(Tables.getcolumn(geom, first(geometrycolumn)))
17+
end
18+
return _is3d(first_geom)
19+
else # assume iterable
20+
first_geom = first(geom)
21+
if GI.trait(first_geom) isa GI.AbstractTrait
22+
return _is3d(first_geom)
23+
else
24+
return false # couldn't figure it out!
25+
end
26+
end
27+
28+
end
829

930
_npoint(x) = _npoint(trait(x), x)
1031
_npoint(::Nothing, xs::AbstractArray) = sum(_npoint, xs)

‎test/primitives.jl

Lines changed: 34 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -56,30 +56,45 @@ poly = GI.Polygon([lr1, lr2])
5656
end
5757
end
5858

59-
@testset "DataFrames" begin
60-
countries_df = DataFrames.DataFrame(countries_table)
61-
GO.DataAPI.metadata!(countries_df, "note metadata", "note metadata value"; style = :note)
62-
GO.DataAPI.metadata!(countries_df, "default metadata", "default metadata value"; style = :default)
63-
centroid_df = GO.apply(GO.centroid, GO.TraitTarget(GI.PolygonTrait(), GI.MultiPolygonTrait()), countries_df; crs = GFT.EPSG(3031));
64-
@test centroid_df isa DataFrames.DataFrame
65-
centroid_geometry = centroid_df.geometry
66-
# Test that the centroids are correct
67-
@test all(centroid_geometry .== GO.centroid.(countries_df.geometry))
68-
@testset "Columns are preserved" begin
69-
for column in Iterators.filter(!=(:geometry), GO.Tables.columnnames(countries_df))
70-
@test all(missing_or_equal.(centroid_df[!, column], countries_df[!, column]))
71-
end
72-
end
73-
@testset "Metadata preservation (or not)" begin
74-
@test DataAPI.metadata(centroid_df, "note metadata") == "note metadata value"
75-
@test !("default metadata" in DataAPI.metadatakeys(centroid_df))
76-
@test DataAPI.metadata(centroid_df, "GEOINTERFACE:geometrycolumns") == (:geometry,)
77-
@test DataAPI.metadata(centroid_df, "GEOINTERFACE:crs") == GFT.EPSG(3031)
59+
@testset "DataFrames" begin
60+
countries_df = DataFrames.DataFrame(countries_table)
61+
GO.DataAPI.metadata!(countries_df, "note metadata", "note metadata value"; style = :note)
62+
GO.DataAPI.metadata!(countries_df, "default metadata", "default metadata value"; style = :default)
63+
centroid_df = GO.apply(GO.centroid, GO.TraitTarget(GI.PolygonTrait(), GI.MultiPolygonTrait()), countries_df; crs = GFT.EPSG(3031));
64+
# Test that the Tables.jl materializer is used
65+
@test centroid_df isa DataFrames.DataFrame
66+
# Test that the centroids are correct
67+
@test all(centroid_df.geometry .== GO.centroid.(countries_df.geometry))
68+
@testset "Columns are preserved" begin
69+
for column in filter(!=(:geometry), GO.Tables.columnnames(countries_df))
70+
@test all(missing_or_equal.(centroid_df[!, column], countries_df[!, column]))
7871
end
7972
end
73+
@testset "Multiple geometry columns in metadata" begin
74+
# set up a dataframe with multiple geometry columns
75+
countries_df2 = deepcopy(countries_df)
76+
countries_df2.centroid = GO.centroid.(countries_df2.geometry)
77+
GI.DataAPI.metadata!(countries_df2, "GEOINTERFACE:geometrycolumns", (:geometry, :centroid); style = :note)
78+
transformed = GO.transform(p -> p .+ 3, countries_df2)
79+
@test GI.DataAPI.metadata(transformed, "GEOINTERFACE:geometrycolumns") == (:geometry, :centroid)
80+
@test GI.DataAPI.metadata(transformed, "GEOINTERFACE:crs") == GI.crs(countries_df2)
81+
# Test that the transformation was actually applied to both geometry columns.
82+
@test all(map(zip(countries_df2.geometry, transformed.geometry)) do (o, n)
83+
GO.equals(GO.transform(p -> p .+ 3, o), n)
84+
end)
85+
@test all(map(zip(countries_df2.centroid, transformed.centroid)) do (o, n)
86+
any(isnan, o) || GO.equals(GO.transform(p -> p .+ 3, o), n)
87+
end)
88+
end
89+
end
8090
end
8191
end
8292
end
93+
@testset "Wrong geometry column kwarg" begin
94+
tab = Tables.dictcolumntable((; geometry = [(1, 2), (3, 4), (5, 6)], other = [1, 2, 3]))
95+
@test_throws "got a Float64" GO.transform(identity, tab; geometrycolumn = 1000.0)
96+
@test_throws "but the table has columns" GO.transform(identity, tab; geometrycolumn = :somethingelse)
97+
end
8398
end
8499
end
85100

‎test/transformations/reproject.jl

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -79,3 +79,16 @@ _xy(p) = GI.x(p), GI.y(p)
7979

8080
GO.reproject(multipolygon4326; target_crs=ProjString("+proj=moll +type=crs"))
8181
end
82+
83+
84+
@testset "Reproject with only target crs" begin
85+
multipolygon4326 = GO.tuples(multipolygon; crs = EPSG(4326))
86+
multipolygon3857 = GO.reproject(multipolygon, EPSG(4326), EPSG(3857))
87+
@test GI.crs(GO.reproject(multipolygon3857; target_crs=EPSG(4326))) == EPSG(4326)
88+
@test GI.crs(GO.reproject(multipolygon3857, EPSG(4326))) == EPSG(4326)
89+
90+
v = [multipolygon3857]
91+
@test GI.crs(only(GO.reproject(v, EPSG(4326)))) == EPSG(4326)
92+
end
93+
94+

0 commit comments

Comments
 (0)