From 08985616c84dacc1affb672d2fbf5245fac794f2 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 8 Aug 2024 11:57:18 -0400 Subject: [PATCH 1/3] Propagate CRS and geometry column from Shapefile.Table via DataAPI Use correct name for geomcols --- Project.toml | 4 +++- src/Shapefile.jl | 2 +- src/table.jl | 35 +++++++++++++++++++++++++++++++++++ test/table.jl | 14 +++++++++++++- 4 files changed, 52 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 2d240c3..e730a5b 100644 --- a/Project.toml +++ b/Project.toml @@ -5,6 +5,7 @@ version = "0.13.0" [deps] DBFTables = "75c7ada1-017a-5fb6-b8c7-2125ff2d6c93" +DataAPI = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a" Extents = "411431e0-e8b7-467b-b5e0-f676ba4f2910" GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f" GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f" @@ -24,6 +25,7 @@ ShapefileZipFileExt = "ZipFile" [compat] DBFTables = "1.2" +DataAPI = "1" Extents = "0.1" GeoFormatTypes = "0.4" GeoInterface = "1.0" @@ -42,8 +44,8 @@ DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" RemoteFiles = "cbe49d4c-5af1-5b60-bb70-0a60aa018e1b" -ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +ZipFile = "a5390f91-8eb1-5f08-bee0-b1d1ffed6cea" [targets] test = ["ArchGDAL", "DataFrames", "Makie", "Plots", "RemoteFiles", "ZipFile", "Test"] diff --git a/src/Shapefile.jl b/src/Shapefile.jl index 97a0abf..7690940 100644 --- a/src/Shapefile.jl +++ b/src/Shapefile.jl @@ -1,6 +1,6 @@ module Shapefile -import GeoFormatTypes, GeoInterface, GeoInterfaceRecipes, DBFTables, Extents, Tables +import GeoFormatTypes, GeoInterface, GeoInterfaceRecipes, DBFTables, Extents, Tables, DataAPI const GI = GeoInterface const GFT = GeoFormatTypes diff --git a/src/table.jl b/src/table.jl index 05bbd64..0b0a7b3 100644 --- a/src/table.jl +++ b/src/table.jl @@ -90,12 +90,47 @@ getdbf(t::Table) = getfield(t, :dbf) Base.length(t::Table) = length(shapes(t)) +# Tables.jl interface Tables.istable(::Type{<:Table}) = true Tables.rows(t::Table) = t Tables.columns(t::Table) = t Tables.rowaccess(::Type{<:Table}) = true Tables.columnaccess(::Type{<:Table}) = true +# DataAPI.jl interface +# Note that DataAPI.jl metadata styles have meaning. +# - `:default` style means that metadata is invalidated when the table is modified. +# - `:note` style means that metadata is preserved when the table is modified. +# Since the CRS and geometry column metadata are properties of the table, and Shapefile.jl +# structures are not meant to be mutable, we use the `:note` style, so that the metadata +# is preserved and propagated correctly. +# NOTE: any `reproject` implementation in e.g. GeometryOps will need to rewrite the CRS metadata +# to a table too, for example. +DataAPI.metadatasupport(::Type{<:Table}) = (; read = true, write = false) +DataAPI.metadatakeys(t::Table) = ("GEOINTERFACE:geometrycolumns", "GEOINTERFACE:crs") + +function DataAPI.metadata(t::Table; style = false) + if style + return Dict("GEOINTERFACE:geometrycolumns" => ((:geometry,), :note), "GEOINTERFACE:crs" => (GeoInterface.crs(t), :note)) + else + return Dict("GEOINTERFACE:geometrycolumns" => (:geometry,), "GEOINTERFACE:crs" => GeoInterface.crs(t)) + end +end +function DataAPI.metadata(t::Table, key::String; style = false) + result = if key == "GEOINTERFACE:geometrycolumns" + (:geometry,) + elseif key == "GEOINTERFACE:crs" + GeoInterface.crs(t) + else + nothing + end + # Check style and return the appropriate type + if style + return (result, :note) + else + return result + end +end """ diff --git a/test/table.jl b/test/table.jl index 1750b55..05d457d 100644 --- a/test/table.jl +++ b/test/table.jl @@ -4,6 +4,7 @@ using RemoteFiles using Plots using Makie import DBFTables +import DataAPI import Tables import DataFrames import GeoInterface @@ -103,14 +104,21 @@ wkt = "GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",637813 (:geometry, :featurecla, :scalerank, :min_zoom), (Union{Missing, Shapefile.Polygon}, Union{String,Missing}, Union{Int,Missing}, Union{Float64,Missing}), ) + # Test DataAPI implementation + @test isempty(setdiff(DataAPI.metadatakeys(ne_land), ("GEOINTERFACE:geometrycolumns", "GEOINTERFACE:crs"))) + @test DataAPI.metadata(ne_land; style = false) isa Dict + @test DataAPI.metadata(ne_land, "GEOINTERFACE:geometrycolumns"; style = false) == (:geometry,) + @test DataAPI.metadata(ne_land, "GEOINTERFACE:crs"; style = false) == GeoInterface.crs(ne_land) for r in ne_land @test Shapefile.shape(r) isa Shapefile.Polygon - @test r.featurecla === "Land" + @test r.featurecla == "Land" end df_land = DataFrames.DataFrame(ne_land) @test size(df_land) == (127, 4) @test names(df_land) == ["geometry", "featurecla", "scalerank", "min_zoom"] df_land.featurecla isa Vector{Union{String,Missing}} + @test DataAPI.metadata(df_land, "GEOINTERFACE:geometrycolumns"; style = false) == (:geometry,) + @test DataAPI.metadata(df_land, "GEOINTERFACE:crs"; style = false) == GeoInterface.crs(df_land) end @testset "ne_coastline" begin @@ -142,6 +150,8 @@ end @test size(df_coastline) == (134, 4) @test names(df_coastline) == ["geometry", "scalerank", "featurecla", "min_zoom"] df_coastline.featurecla isa Vector{Union{String,Missing}} + @test DataAPI.metadata(df_coastline, "GEOINTERFACE:geometrycolumns"; style = false) == (:geometry,) + @test DataAPI.metadata(df_coastline, "GEOINTERFACE:crs"; style = false) == GeoInterface.crs(df_coastline) end @testset "ne_cities" begin @@ -187,6 +197,8 @@ end @test size(df_cities) == (243, 39) @test names(df_cities) == string.(colnames) df_cities.featurecla isa Vector{Union{String,Missing}} + @test DataAPI.metadata(df_cities, "GEOINTERFACE:geometrycolumns"; style = false) == (:geometry,) + @test DataAPI.metadata(df_cities, "GEOINTERFACE:crs"; style = false) == GeoInterface.crs(df_cities) end # no need to use shx in Shapefile.Tables since we read the shapes into a Vector and can thus index them From b1398ad458e182249544823758462ef355624d70 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 8 Aug 2024 12:04:14 -0400 Subject: [PATCH 2/3] Fix tests for now since GeoInterface doesn't subscribe to DataAPI yet --- test/table.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/table.jl b/test/table.jl index 05d457d..1ac941f 100644 --- a/test/table.jl +++ b/test/table.jl @@ -118,7 +118,7 @@ wkt = "GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",637813 @test names(df_land) == ["geometry", "featurecla", "scalerank", "min_zoom"] df_land.featurecla isa Vector{Union{String,Missing}} @test DataAPI.metadata(df_land, "GEOINTERFACE:geometrycolumns"; style = false) == (:geometry,) - @test DataAPI.metadata(df_land, "GEOINTERFACE:crs"; style = false) == GeoInterface.crs(df_land) + @test DataAPI.metadata(df_land, "GEOINTERFACE:crs"; style = false) == GeoInterface.crs(ne_land) end @testset "ne_coastline" begin @@ -151,7 +151,7 @@ end @test names(df_coastline) == ["geometry", "scalerank", "featurecla", "min_zoom"] df_coastline.featurecla isa Vector{Union{String,Missing}} @test DataAPI.metadata(df_coastline, "GEOINTERFACE:geometrycolumns"; style = false) == (:geometry,) - @test DataAPI.metadata(df_coastline, "GEOINTERFACE:crs"; style = false) == GeoInterface.crs(df_coastline) + @test DataAPI.metadata(df_coastline, "GEOINTERFACE:crs"; style = false) == GeoInterface.crs(ne_coastline) end @testset "ne_cities" begin @@ -198,7 +198,7 @@ end @test names(df_cities) == string.(colnames) df_cities.featurecla isa Vector{Union{String,Missing}} @test DataAPI.metadata(df_cities, "GEOINTERFACE:geometrycolumns"; style = false) == (:geometry,) - @test DataAPI.metadata(df_cities, "GEOINTERFACE:crs"; style = false) == GeoInterface.crs(df_cities) + @test DataAPI.metadata(df_cities, "GEOINTERFACE:crs"; style = false) == GeoInterface.crs(ne_cities) end # no need to use shx in Shapefile.Tables since we read the shapes into a Vector and can thus index them From 060843fd7a626e60fef630a24c74921ca23b5698 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 8 Aug 2024 12:04:34 -0400 Subject: [PATCH 3/3] Try to get CRS from DataAPI metadata in writer --- src/writer.jl | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/writer.jl b/src/writer.jl index 1966a06..54df654 100644 --- a/src/writer.jl +++ b/src/writer.jl @@ -34,12 +34,19 @@ struct Writer all(x -> GI.isgeometry(x) || ismissing(x), geoms) || error("Not all geoms satisfy `GeoInterface.isgeometry`.") + # Try to get the CRS from the table if it exists + writable_crs = if isnothing(crs) && DataAPI.metadatasupport(feats).read && "GEOINTERFACE:crs" in DataAPI.metadatakeys(feats) + DataAPI.metadata(feats, "GEOINTERFACE:crs"; style = false) + else + crs + end + ngeoms = sum(1 for _ in geoms) nfeats = sum(1 for _ in Tables.rows(feats)) ngeoms == nfeats || error("Number of geoms does not match number of features. Found: $ngeoms ≠ $nfeats.") - new(geoms, feats, crs) + new(geoms, feats, writable_crs) end end