Propagate CRS and geometry column from Shapefile.Table via DataAPI

Project.toml

@@ -5,6 +5,7 @@ version = "0.13.1"
 
 [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"

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
 

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
 
 
 """

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
 

test/table.jl

@@ -4,6 +4,7 @@ using RemoteFiles
 using Plots
 using Makie
 import DBFTables
+import DataAPI
 import Tables
 import DataFrames
 import GeoInterface
@@ -103,6 +104,11 @@ 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"
@@ -111,6 +117,8 @@ wkt = "GEOGCS[\"GCS_WGS_1984\",DATUM[\"D_WGS_1984\",SPHEROID[\"WGS_1984\",637813
  @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(ne_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(ne_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(ne_cities)
 end
 
 # no need to use shx in Shapefile.Tables since we read the shapes into a Vector and can thus index them