Support GeoInterface geometry column(s). (#41)

* Support GeoInterface geometry column(s).

* Update dependend on recent GeoInterface.

* Always try to write :geometry as column name.

* Catch geom_columns being nothing.

* Up GI compat.

.gitignore

@@ -9,3 +9,5 @@ docs/src/*.png
 *.sbn
 *.sbx
 *.sha256
+.DS_Store
+Manifest.toml

Project.toml

@@ -1,19 +1,21 @@
 name = "GeoDataFrames"
 uuid = "62cb38b5-d8d2-4862-a48e-6a340996859f"
 authors = ["Maarten Pronk <[email protected]> and contributors"]
-version = "0.2.4"
+version = "0.3.0"
 
 [deps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 GeoFormatTypes = "68eda718-8dee-11e9-39e7-89f7f65f511f"
+GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
 [compat]
-ArchGDAL = "0.8, 0.9"
+ArchGDAL = "0.9"
 DataFrames = "1.0"
 GeoFormatTypes = "0.3, 0.4"
+GeoInterface = "1.0.1"
 Tables = "1"
 julia = "1.6"
 

src/GeoDataFrames.jl

@@ -4,8 +4,10 @@ import ArchGDAL as AG
 using DataFrames
 using Tables
 import GeoFormatTypes as GFT
+import GeoInterface
 
 include("exports.jl")
 include("io.jl")
+include("utils.jl")
 
 end  # module

src/io.jl

@@ -10,6 +10,17 @@ const drivermapping = Dict(
     ".nc" => "netCDF",
 )
 
+const lookup_type = Dict{DataType,AG.OGRwkbGeometryType}(
+    AG.GeoInterface.PointTrait => AG.wkbPoint,
+    AG.GeoInterface.MultiPointTrait => AG.wkbMultiPoint,
+    AG.GeoInterface.LineStringTrait => AG.wkbLineString,
+    AG.GeoInterface.LinearRingTrait => AG.wkbMultiLineString,
+    AG.GeoInterface.MultiLineStringTrait => AG.wkbMultiLineString,
+    AG.GeoInterface.PolygonTrait => AG.wkbPolygon,
+    AG.GeoInterface.MultiPolygonTrait => AG.wkbMultiPolygon,
+)
+
+
 """
     read(fn::AbstractString; kwargs...)
     read(fn::AbstractString, layer::Union{Integer,AbstractString}; kwargs...)
@@ -41,21 +52,37 @@ function read(ds, layer)
         end
         return DataFrame(table)
     end
-    "" in names(df) && rename!(df, Dict(Symbol("") => :geom,))  # needed for now
+    "" in names(df) && rename!(df, Dict(Symbol("") => :geometry,))  # needed for now
     return df
 end
 
 """
-    write(fn::AbstractString, table; layer_name="data", geom_column=:geom, crs::Union{GFT.GeoFormat,Nothing}=nothing, driver::Union{Nothing,AbstractString}=nothing)
+    write(fn::AbstractString, table; layer_name="data", geom_column=:geometry, crs::Union{GFT.GeoFormat,Nothing}=nothing, driver::Union{Nothing,AbstractString}=nothing, options::Vector{AbstractString}=[], geom_columns::Set{Symbol}=(:geometry))
 
 Write the provided `table` to `fn`. The `geom_column` is expected to hold ArchGDAL geometries.
 """
-function write(fn::AbstractString, table; layer_name::AbstractString="data", geom_column::Symbol=:geom, crs::Union{GFT.GeoFormat,Nothing}=nothing, driver::Union{Nothing,AbstractString}=nothing)
+function write(fn::AbstractString, table; layer_name::AbstractString="data", crs::Union{GFT.GeoFormat,Nothing}=nothing, driver::Union{Nothing,AbstractString}=nothing, options::Dict{String,String}=Dict{String,String}(), geom_columns=GeoInterface.geometrycolumns(table), kwargs...)
     rows = Tables.rows(table)
     sch = Tables.schema(rows)
 
-    # Geom type can't be inferred from here
-    geom_type = AG.getgeomtype(rows[1][geom_column])
+    # Determine geometry columns
+    isnothing(geom_columns) && error("Please set `geom_columns` or define `GeoInterface.geometrycolumns` for $(typeof(table))")
+    if :geom_column in keys(kwargs)  # backwards compatible
+        geom_columns = (kwargs[:geom_column],)
+    end
+
+    geom_types = []
+    for geom_column in geom_columns
+        trait = AG.GeoInterface.geomtrait(getproperty(first(rows), geom_column))
+        geom_type = get(lookup_type, typeof(trait), nothing)
+        isnothing(geom_type) && throw(ArgumentError("Can't convert $trait of column $geom_column to ArchGDAL yet."))
+        push!(geom_types, geom_type)
+    end
+
+    # Set geometry name in options
+    if !("geometry_name" in keys(options))
+        options["geometry_name"] = "geometry"
+    end
 
     # Find driver
     _, extension = splitext(fn)
@@ -70,14 +97,15 @@ function write(fn::AbstractString, table; layer_name::AbstractString="data", geo
     # Figure out attributes
     fields = Vector{Tuple{Symbol,DataType}}()
     for (name, type) in zip(sch.names, sch.types)
-        if type != AG.IGeometry && name != geom_column
+        if !(name in geom_columns)
+            AG.GeoInterface.isgeometry(type) && error("Did you mean to use the `geom_columns` argument to specify $name is a geometry?")
             types = Base.uniontypes(type)
             if length(types) == 1
                 push!(fields, (Symbol(name), type))
             elseif length(types) == 2 && Missing in types
                 push!(fields, (Symbol(name), types[2]))
             else
-                error("Can't convert to GDAL type from $type")
+                error("Can't convert to GDAL type from $type. Please file an issue.")
             end
         end
     end
@@ -89,9 +117,15 @@ function write(fn::AbstractString, table; layer_name::AbstractString="data", geo
             crs !== nothing && AG.importCRS!(spatialref, crs)
             AG.createlayer(
                 name=layer_name,
-                geom=geom_type,
-                spatialref=spatialref
+                geom=first(geom_types),  # how to set the name though?
+                spatialref=spatialref,
+                options=stringlist(options)
             ) do layer
+                for (i, (geom_column, geom_type)) in enumerate(zip(geom_columns, geom_types))
+                    if i > 1
+                        AG.writegeomdefn!(layer, string(geom_column), geom_type)
+                    end
+                end
                 for (name, type) in fields
                     AG.createfielddefn(String(name), convert(AG.OGRFieldType, type)) do fd
                         AG.setsubtype!(fd, convert(AG.OGRFieldSubType, type))
@@ -100,7 +134,9 @@ function write(fn::AbstractString, table; layer_name::AbstractString="data", geo
                 end
                 for row in rows
                     AG.createfeature(layer) do feature
-                        AG.setgeom!(feature, getproperty(row, geom_column))
+                        for (i, (geom_column)) in enumerate(geom_columns)
+                            AG.setgeom!(feature, i - 1, convert(AG.IGeometry, getproperty(row, geom_column)))
+                        end
                         for (name, _) in fields
                             field = getproperty(row, name)
                             if !ismissing(field)
@@ -111,7 +147,7 @@ function write(fn::AbstractString, table; layer_name::AbstractString="data", geo
                         end
                     end
                 end
-                AG.copy(layer, dataset=ds, name=layer_name)
+                AG.copy(layer, dataset=ds, name=layer_name, options=stringlist(options))
             end
         end
     end

src/utils.jl

@@ -0,0 +1,7 @@
+function stringlist(dict::Dict{String,String})
+    sv = Vector{String}()
+    for (k, v) in pairs(dict)
+        push!(sv, uppercase(string(k)) * "=" * string(v))
+    end
+    return sv
+end

test/runtests.jl

@@ -1,10 +1,11 @@
+import GeoDataFrames as GDF
 using DataFrames
 using Dates
 using Pkg.PlatformEngines
 using Test
 import ArchGDAL as AG
-import GeoDataFrames as GDF
 import GeoFormatTypes as GFT
+import GeoInterface as GI
 
 # Use ArchGDAL datasets to test with
 const testdatadir = joinpath(@__DIR__, "data")
@@ -64,10 +65,10 @@ end
 
     @testset "Read self written file" begin
         # Save table with a few random points
-        table = DataFrame(geom=AG.createpoint.(coords), name="test")
+        table = DataFrame(geometry=AG.createpoint.(coords), name="test")
         GDF.write(joinpath(testdatadir, "test_points.shp"), table)
-        GDF.write(joinpath(testdatadir, "test_points.gpkg"), table, layer_name="test_points")
-        GDF.write(joinpath(testdatadir, "test_points.geojson"), table, layer_name="test_points")
+        GDF.write(joinpath(testdatadir, "test_points.gpkg"), table; layer_name="test_points")
+        GDF.write(joinpath(testdatadir, "test_points.geojson"), table; layer_name="test_points")
 
         ntable = GDF.read(joinpath(testdatadir, "test_points.shp"))
         @test nrow(ntable) == 10
@@ -78,20 +79,19 @@ end
     end
 
     @testset "Write shapefile" begin
-
         t = GDF.read(fn)
 
         # Save table from reading
-        GDF.write(joinpath(testdatadir, "test_read.shp"), t, layer_name="test_coastline")
-        GDF.write(joinpath(testdatadir, "test_read.gpkg"), t, layer_name="test_coastline")
-        GDF.write(joinpath(testdatadir, "test_read.geojson"), t, layer_name="test_coastline")
+        GDF.write(joinpath(testdatadir, "test_read.shp"), t; layer_name="test_coastline")
+        GDF.write(joinpath(testdatadir, "test_read.gpkg"), t; layer_name="test_coastline")
+        GDF.write(joinpath(testdatadir, "test_read.geojson"), t; layer_name="test_coastline")
 
     end
 
     @testset "Write shapefile with non-GDAL types" begin
         coords = collect(zip(rand(Float32, 2), rand(Float32, 2)))
         t = DataFrame(
-            geom=AG.createpoint.(coords),
+            geometry=AG.createpoint.(coords),
             name=["test", "test2"],
             flag=UInt8[typemin(UInt8), typemax(UInt8)],
             ex1=Int16[typemin(Int8), typemax(Int8)],
@@ -107,7 +107,7 @@ end
         GDF.write(joinpath(testdatadir, "test_exotic.gpkg"), t)
         GDF.write(joinpath(testdatadir, "test_exotic.geojson"), t)
         tt = GDF.read(joinpath(testdatadir, "test_exotic.gpkg"))
-        @test AG.getx.(tt.geom, 0) == AG.getx.(t.geom, 0)
+        @test AG.getx.(tt.geometry, 0) == AG.getx.(t.geometry, 0)
         @test tt.flag == t.flag
         @test tt.ex1 == t.ex1
         @test tt.ex2 == t.ex2
@@ -126,21 +126,43 @@ end
     end
 
     @testset "Spatial operations" begin
-        table = DataFrame(geom=AG.createpoint.(coords), name="test")
+        table = DataFrame(geometry=AG.createpoint.(coords), name="test")
 
         # Buffer to also write polygons
-        table.geom = AG.buffer(table.geom, 10)
+        table.geometry = AG.buffer(table.geometry, 10)
         GDF.write(joinpath(testdatadir, "test_polygons.shp"), table)
         GDF.write(joinpath(testdatadir, "test_polygons.gpkg"), table)
         GDF.write(joinpath(testdatadir, "test_polygons.geojson"), table)
     end
 
     @testset "Reproject" begin
-        table = DataFrame(geom=AG.createpoint.([[0, 0, 0]]), name="test")
-        AG.reproject(table.geom, GFT.EPSG(4326), GFT.EPSG(28992))
-        @test GDF.AG.getpoint(table.geom[1], 0)[1] ≈ -587791.596556932
-        GDF.write(joinpath(testdatadir, "test_reprojection.gpkg"), table, crs=GFT.EPSG(28992))
+        table = DataFrame(geometry=AG.createpoint.([[0, 0, 0]]), name="test")
+        AG.reproject(table.geometry, GFT.EPSG(4326), GFT.EPSG(28992))
+        @test GDF.AG.getpoint(table.geometry[1], 0)[1] ≈ -587791.596556932
+        GDF.write(joinpath(testdatadir, "test_reprojection.gpkg"), table; crs=GFT.EPSG(28992))
     end
 
+    @testset "Kwargs" begin
+        table = DataFrame(foo=AG.createpoint.([[0, 0, 0]]), name="test")
+        GDF.write(joinpath(testdatadir, "test_options1.gpkg"), table; geom_column=:foo)
+        GDF.write(joinpath(testdatadir, "test_options2.gpkg"), table; geom_columns=Set((:foo,)))
+
+        table = DataFrame(foo=AG.createpoint.([[0, 0, 0]]), bar=AG.createpoint.([[0, 0, 0]]), name="test")
+        @test_throws Exception GDF.write(joinpath(testdatadir, "test_options3.gpkg"), table; geom_column=:foo)
+        GDF.write(joinpath(testdatadir, "test_options3.gpkg"), table; geom_columns=Set((:foo, :bar)))
+
+        table = DataFrame(foo=AG.createpoint.([[0, 0, 0]]), name="test")
+        GDF.write(joinpath(testdatadir, "test_options4.gpkg"), table; options=Dict(
+                "GEOMETRY_NAME" => "bar", "DESCRIPTION" => "Written by GeoDataFrames.jl"), geom_column=:foo)
+
+    end
+
+    @testset "GeoInterface" begin
+        tfn = joinpath(testdatadir, "test_geointerface.gpkg")
+        table = [(; geom=AG.createpoint(1.0, 2.0), name="test")]
+        @test_throws Exception GDF.write(tfn, table)
+        GI.geometrycolumns(::Vector{<:NamedTuple}) = (:geom,)
+        @test isfile(GDF.write(tfn, table))
+    end
 
 end