Add native driver extensions

Project.toml

@@ -14,14 +14,30 @@ Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
+[weakdeps]
+FlatGeobuf = "d985ece1-97de-4d33-914c-38fb84042e15"
+GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9"
+GeoParquet = "e99870d8-ce00-4fdd-aeee-e09192881159"
+Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
+
+[extensions]
+GeoDataFramesFlatGeobufExt = "FlatGeobuf"
+GeoDataFramesGeoJSONExt = "GeoJSON"
+GeoDataFramesGeoParquetExt = "GeoParquet"
+GeoDataFramesShapefileExt = "Shapefile"
+
 [compat]
 ArchGDAL = "0.10"
 DataAPI = "1.13"
 DataFrames = "1.4"
 Extents = "0.1"
+FlatGeobuf = "0.1.2"
 GeoFormatTypes = "0.3, 0.4"
 GeoInterface = "1.0.1"
+GeoJSON = "0.8.1"
+GeoParquet = "0.2.1"
 Reexport = "1.2"
+Shapefile = "0.13.1"
 Tables = "1"
 julia = "1.9"
 
@@ -30,4 +46,4 @@ Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Dates"]
+test = ["Test", "Dates", "Shapefile", "GeoJSON", "GeoParquet", "FlatGeobuf"]

ext/GeoDataFramesFlatGeobufExt.jl

@@ -0,0 +1,16 @@
+module GeoDataFramesFlatGeobufExt
+
+using FlatGeobuf
+using GeoDataFrames: FlatGeobufDriver, ArchGDALDriver, GeoDataFrames
+
+function GeoDataFrames.read(::FlatGeobufDriver, fname::AbstractString; kwargs...)
+    df = GeoDataFrames.DataFrame(FlatGeobuf.read_file(fname; kwargs...))
+    GeoDataFrames.rename!(df, :geom => :geometry)
+    return df
+end
+
+function GeoDataFrames.write(::FlatGeobufDriver, fname::AbstractString, data; kwargs...)
+    # No write support yet
+    GeoDataFrames.write(ArchGDALDriver(), fname, data; kwargs...)
+end
+end

ext/GeoDataFramesGeoJSONExt.jl

@@ -0,0 +1,14 @@
+module GeoDataFramesGeoJSONExt
+
+using GeoDataFrames: GeoJSONDriver, GeoDataFrames
+using GeoJSON
+
+function GeoDataFrames.read(::GeoJSONDriver, fname::AbstractString; kwargs...)
+    GeoDataFrames.DataFrame(GeoJSON.read(fname; kwargs...))
+end
+
+function GeoDataFrames.write(::GeoJSONDriver, fname::AbstractString, data; kwargs...)
+    GeoJSON.write(fname, data; kwargs...)
+end
+
+end

ext/GeoDataFramesGeoParquetExt.jl

@@ -0,0 +1,19 @@
+module GeoDataFramesGeoParquetExt
+
+using GeoDataFrames: GeoParquetDriver, GeoDataFrames
+import GeoInterface as GI
+using GeoParquet
+
+function GeoDataFrames.read(::GeoParquetDriver, fname::AbstractString; kwargs...)
+    df = GeoParquet.read(fname; kwargs...)
+    crs = GeoDataFrames.metadata(df, "GEOINTERFACE:crs")
+    ncrs =  GeoDataFrames.GFT.ProjJSON(GeoParquet.JSON3.write(crs.val))
+    GeoDataFrames.metadata!(df, "GEOINTERFACE:crs", ncrs; style = :note)
+    df
+end
+
+function GeoDataFrames.write(::GeoParquetDriver, fname::AbstractString, data; kwargs...)
+    GeoParquet.write(fname, data, GI.geometrycolumns(data); kwargs...)
+end
+
+end

ext/GeoDataFramesShapefileExt.jl

@@ -0,0 +1,14 @@
+module GeoDataFramesShapefileExt
+
+using GeoDataFrames: ShapefileDriver, GeoDataFrames
+using Shapefile
+
+function GeoDataFrames.read(::ShapefileDriver, fname::AbstractString; kwargs...)
+    GeoDataFrames.DataFrame(Shapefile.Table(fname; kwargs...); copycols = false)
+end
+
+function GeoDataFrames.write(::ShapefileDriver, fname::AbstractString, data; kwargs...)
+    Shapefile.write(fname, data; kwargs...)
+end
+
+end

src/GeoDataFrames.jl

@@ -4,12 +4,16 @@ import ArchGDAL as AG
 using DataFrames
 using Tables
 import GeoFormatTypes as GFT
-import GeoInterface
+import GeoInterface as GI
 using DataAPI
 using Reexport
 
 include("exports.jl")
+include("drivers.jl")
 include("io.jl")
 include("utils.jl")
 
+function load end
+function save end
+
 end  # module

src/drivers.jl

@@ -0,0 +1,33 @@
+abstract type AbstractDriver end
+
+struct GeoJSONDriver <: AbstractDriver end
+struct ShapefileDriver <: AbstractDriver end
+struct GeoParquetDriver <: AbstractDriver end
+struct FlatGeobufDriver <: AbstractDriver end
+struct ArchGDALDriver <: AbstractDriver end
+
+function driver(ext::AbstractString)
+    if ext in (".json", ".geojson")
+        return GeoJSONDriver()
+    elseif ext == ".shp"
+        return ShapefileDriver()
+    elseif ext in (".parquet", ".pq")
+        return GeoParquetDriver()
+    elseif ext == ".fgb"
+        return FlatGeobufDriver()
+    else
+        return ArchGDALDriver()
+    end
+end
+
+package(::GeoJSONDriver) = :GeoJSON
+package(::ShapefileDriver) = :Shapefile
+package(::GeoParquetDriver) = :GeoParquet
+package(::FlatGeobufDriver) = :FlatGeobuf
+package(::ArchGDALDriver) = :ArchGDAL
+
+uuid(::GeoJSONDriver) = "61d90e0f-e114-555e-ac52-39dfb47a3ef9"
+uuid(::ShapefileDriver) = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4"
+uuid(::GeoParquetDriver) = "e99870d8-ce00-4fdd-aeee-e09192881159"
+uuid(::FlatGeobufDriver) = "d985ece1-97de-4d33-914c-38fb84042e15"
+uuid(::ArchGDALDriver) = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"

src/io.jl

@@ -22,24 +22,24 @@ function find_driver(fn::AbstractString)
 end
 
 const lookup_type = Dict{Tuple{DataType, Int}, AG.OGRwkbGeometryType}(
-    (AG.GeoInterface.PointTrait, 2) => AG.wkbPoint,
-    (AG.GeoInterface.PointTrait, 3) => AG.wkbPoint25D,
-    (AG.GeoInterface.PointTrait, 4) => AG.wkbPointZM,
-    (AG.GeoInterface.MultiPointTrait, 2) => AG.wkbMultiPoint,
-    (AG.GeoInterface.MultiPointTrait, 3) => AG.wkbMultiPoint25D,
-    (AG.GeoInterface.MultiPointTrait, 4) => AG.wkbMultiPointZM,
-    (AG.GeoInterface.LineStringTrait, 2) => AG.wkbLineString,
-    (AG.GeoInterface.LineStringTrait, 3) => AG.wkbLineString25D,
-    (AG.GeoInterface.LineStringTrait, 4) => AG.wkbLineStringZM,
-    (AG.GeoInterface.MultiLineStringTrait, 2) => AG.wkbMultiLineString,
-    (AG.GeoInterface.MultiLineStringTrait, 3) => AG.wkbMultiLineString25D,
-    (AG.GeoInterface.MultiLineStringTrait, 4) => AG.wkbMultiLineStringZM,
-    (AG.GeoInterface.PolygonTrait, 2) => AG.wkbPolygon,
-    (AG.GeoInterface.PolygonTrait, 3) => AG.wkbPolygon25D,
-    (AG.GeoInterface.PolygonTrait, 4) => AG.wkbPolygonZM,
-    (AG.GeoInterface.MultiPolygonTrait, 2) => AG.wkbMultiPolygon,
-    (AG.GeoInterface.MultiPolygonTrait, 3) => AG.wkbMultiPolygon25D,
-    (AG.GeoInterface.MultiPolygonTrait, 4) => AG.wkbMultiPolygonZM,
+    (GI.PointTrait, 2) => AG.wkbPoint,
+    (GI.PointTrait, 3) => AG.wkbPoint25D,
+    (GI.PointTrait, 4) => AG.wkbPointZM,
+    (GI.MultiPointTrait, 2) => AG.wkbMultiPoint,
+    (GI.MultiPointTrait, 3) => AG.wkbMultiPoint25D,
+    (GI.MultiPointTrait, 4) => AG.wkbMultiPointZM,
+    (GI.LineStringTrait, 2) => AG.wkbLineString,
+    (GI.LineStringTrait, 3) => AG.wkbLineString25D,
+    (GI.LineStringTrait, 4) => AG.wkbLineStringZM,
+    (GI.MultiLineStringTrait, 2) => AG.wkbMultiLineString,
+    (GI.MultiLineStringTrait, 3) => AG.wkbMultiLineString25D,
+    (GI.MultiLineStringTrait, 4) => AG.wkbMultiLineStringZM,
+    (GI.PolygonTrait, 2) => AG.wkbPolygon,
+    (GI.PolygonTrait, 3) => AG.wkbPolygon25D,
+    (GI.PolygonTrait, 4) => AG.wkbPolygonZM,
+    (GI.MultiPolygonTrait, 2) => AG.wkbMultiPolygon,
+    (GI.MultiPolygonTrait, 3) => AG.wkbMultiPolygon25D,
+    (GI.MultiPolygonTrait, 4) => AG.wkbMultiPolygonZM,
 )
 
 """
@@ -49,35 +49,41 @@ const lookup_type = Dict{Tuple{DataType, Int}, AG.OGRwkbGeometryType}(
 Read a file into a DataFrame. Any kwargs are passed onto ArchGDAL [here](https://yeesian.com/ArchGDAL.jl/stable/reference/#ArchGDAL.read-Tuple{AbstractString}).
 By default you only get the first layer, unless you specify either the index (0 based) or name (string) of the layer.
 """
-function read(fn::AbstractString; kwargs...)
+function read(fn; kwargs...)
+    ext = last(splitext(fn))
+    dr = driver(ext)
+    read(dr, fn; kwargs...)
+end
+
+function read(driver::AbstractDriver, fn::AbstractString; kwargs...)
+    @info "Using GDAL for reading, import $(package(driver)) for a native driver."
+    read(ArchGDALDriver(), fn; kwargs...)
+end
+
+function read(driver::ArchGDALDriver, fn::AbstractString; layer=nothing, kwargs...)
     startswith(fn, "/vsi") ||
         occursin(":", fn) ||
         isfile(fn) ||
         isdir(fn) ||
         error("File not found.")
+
     t = AG.read(fn; kwargs...) do ds
         ds.ptr == C_NULL && error("Unable to open $fn.")
-        if AG.nlayer(ds) > 1
-            @warn "This file has multiple layers, you only get the first layer by default now."
+        if AG.nlayer(ds) > 1 && isnothing(layer)
+            @warn "This file has multiple layers, defaulting to first layer."
         end
-        return read(ds, 0)
+        return read(driver, ds, isnothing(layer) ? 0 : layer)
     end
     return t
 end
 
-function read(fn::AbstractString, layer::Union{Integer, AbstractString}; kwargs...)
-    startswith(fn, "/vsi") ||
-        occursin(":", fn) ||
-        isfile(fn) ||
-        isdir(fn) ||
-        error("File not found: $fn")
-    t = AG.read(fn; kwargs...) do ds
-        return read(ds, layer)
-    end
-    return t
-end
+@deprecate read(fn::AbstractString, layer::Union{AbstractString, Integer}; kwargs...) read(
+    fn;
+    layer,
+    kwargs...,
+)
 
-function read(ds, layer)
+function read(::ArchGDALDriver, ds, layer)
     df, gnames, sr = AG.getlayer(ds, layer) do table
         if table.ptr == C_NULL
             throw(
@@ -93,6 +99,9 @@ function read(ds, layer)
     if "" in names(df)
         rename!(df, Symbol("") => :geometry)
         replace!(gnames, Symbol("") => :geometry)
+    elseif "geom" in names(df)
+        rename!(df, Symbol("geom") => :geometry)
+        replace!(gnames, Symbol("geom") => :geometry)
     end
     crs = sr.ptr == C_NULL ? nothing : GFT.WellKnownText(GFT.CRS(), AG.toWKT(sr))
     geometrycolumns = Tuple(gnames)
@@ -110,7 +119,18 @@ end
 
 Write the provided `table` to `fn`. The `geom_column` is expected to hold ArchGDAL geometries.
 """
+function write(fn::AbstractString, table; kwargs...)
+    ext = last(splitext(fn))
+    write(driver(ext), fn, table; kwargs...)
+end
+
+function write(driver::AbstractDriver, fn::AbstractString, table; kwargs...)
+    @info "Using GDAL for writing, import $(package(driver)) for a native driver."
+    write(ArchGDALDriver(), fn, table; kwargs...)
+end
+
 function write(
+    ::ArchGDALDriver,
     fn::AbstractString,
     table;
     layer_name::AbstractString = "data",
@@ -126,16 +146,17 @@ function write(
 
     # Determine geometry columns
     isnothing(geom_columns) && error(
-        "Please set `geom_columns` kw or define `GeoInterface.geometrycolumns` for $(typeof(table))",
+        "Please set `geom_columns` kw or define `GI.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))
-        ndim = AG.GeoInterface.ncoord(getproperty(first(rows), geom_column))
+        geometry = getproperty(first(rows), geom_column)
+        trait = GI.geomtrait(geometry)
+        ndim = GI.ncoord(geometry)
         geom_type = get(lookup_type, (typeof(trait), ndim), nothing)
         isnothing(geom_type) && throw(
             ArgumentError(
@@ -161,7 +182,7 @@ function write(
     fields = Vector{Tuple{Symbol, DataType}}()
     for (name, type) in zip(sch.names, sch.types)
         if !(name in geom_columns)
-            AG.GeoInterface.isgeometry(type) &&
+            GI.isgeometry(type) &&
                 @warn "Writing $name as a non-spatial column, use the `geom_columns` argument to write as a geometry."
             nmtype = nonmissingtype(type)
             if !hasmethod(convert, (Type{AG.OGRFieldType}, Type{nmtype}))
@@ -259,21 +280,21 @@ end
 
 # This should be upstreamed to ArchGDAL
 const lookup_method = Dict{DataType, Function}(
-    GeoInterface.PointTrait => AG.unsafe_createpoint,
-    GeoInterface.MultiPointTrait => AG.unsafe_createmultipoint,
-    GeoInterface.LineStringTrait => AG.unsafe_createlinestring,
-    GeoInterface.LinearRingTrait => AG.unsafe_createlinearring,
-    GeoInterface.MultiLineStringTrait => AG.unsafe_createmultilinestring,
-    GeoInterface.PolygonTrait => AG.unsafe_createpolygon,
-    GeoInterface.MultiPolygonTrait => AG.unsafe_createmultipolygon,
+    GI.PointTrait => AG.unsafe_createpoint,
+    GI.MultiPointTrait => AG.unsafe_createmultipoint,
+    GI.LineStringTrait => AG.unsafe_createlinestring,
+    GI.LinearRingTrait => AG.unsafe_createlinearring,
+    GI.MultiLineStringTrait => AG.unsafe_createmultilinestring,
+    GI.PolygonTrait => AG.unsafe_createpolygon,
+    GI.MultiPolygonTrait => AG.unsafe_createmultipolygon,
 )
 
 function _convert(::Type{T}, geom) where {T <: AG.Geometry}
-    f = get(lookup_method, typeof(GeoInterface.geomtrait(geom)), nothing)
+    f = get(lookup_method, typeof(GI.geomtrait(geom)), nothing)
     isnothing(f) && error(
         "Cannot convert an object of $(typeof(geom)) with the $(typeof(type)) trait (yet). Please report an issue.",
     )
-    return f(GeoInterface.coordinates(geom))
+    return f(GI.coordinates(geom))
 end
 
 function _convert(::Type{T}, geom::AG.IGeometry) where {T <: AG.Geometry}