diff --git a/Project.toml b/Project.toml index 03f15b0..86abc61 100644 --- a/Project.toml +++ b/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"] diff --git a/ext/GeoDataFramesFlatGeobufExt.jl b/ext/GeoDataFramesFlatGeobufExt.jl new file mode 100644 index 0000000..8ff356c --- /dev/null +++ b/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 diff --git a/ext/GeoDataFramesGeoJSONExt.jl b/ext/GeoDataFramesGeoJSONExt.jl new file mode 100644 index 0000000..8db09eb --- /dev/null +++ b/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 diff --git a/ext/GeoDataFramesGeoParquetExt.jl b/ext/GeoDataFramesGeoParquetExt.jl new file mode 100644 index 0000000..278e2b6 --- /dev/null +++ b/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 diff --git a/ext/GeoDataFramesShapefileExt.jl b/ext/GeoDataFramesShapefileExt.jl new file mode 100644 index 0000000..f072024 --- /dev/null +++ b/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 diff --git a/src/GeoDataFrames.jl b/src/GeoDataFrames.jl index b5fccdd..c9617a0 100644 --- a/src/GeoDataFrames.jl +++ b/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 diff --git a/src/drivers.jl b/src/drivers.jl new file mode 100644 index 0000000..53b3bf6 --- /dev/null +++ b/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" diff --git a/src/io.jl b/src/io.jl index 9883646..1530462 100644 --- a/src/io.jl +++ b/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,7 +146,7 @@ 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],) @@ -134,8 +154,9 @@ function write( 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} diff --git a/src/utils.jl b/src/utils.jl index 7ad7e5d..6195c76 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -7,8 +7,8 @@ function stringlist(dict::Dict{String, String}) end function getgeometrycolumns(table) - if GeoInterface.isfeaturecollection(table) - return GeoInterface.geometrycolumns(table) + if GI.isfeaturecollection(table) + return GI.geometrycolumns(table) elseif first(DataAPI.metadatasupport(typeof(table))) gc = DataAPI.metadata(table, "GEOINTERFACE:geometrycolumns", nothing) if isnothing(gc) # fall back to searching for "geometrycolumns" as a string @@ -21,8 +21,8 @@ function getgeometrycolumns(table) end function getcrs(table) - if GeoInterface.isfeaturecollection(table) - return GeoInterface.crs(table) + if GI.isfeaturecollection(table) + return GI.crs(table) end crs = geomcrs(table) if !isnothing(crs) @@ -42,7 +42,7 @@ function geomcrs(table) rows = Tables.rows(table) geom_column = first(getgeometrycolumns(table)) if hasproperty(first(rows), geom_column) - return GeoInterface.crs(getproperty(first(rows), geom_column)) + return GI.crs(getproperty(first(rows), geom_column)) else return nothing end @@ -52,7 +52,7 @@ end # These are the basic metadata definitions from which all else follows. -function GeoInterface.crs(table::DataFrame) +function GI.crs(table::DataFrame) crs = DataAPI.metadata(table, "GEOINTERFACE:crs", nothing) if isnothing(crs) # fall back to searching for "crs" as a string crs = DataAPI.metadata(table, "crs", nothing) @@ -60,7 +60,7 @@ function GeoInterface.crs(table::DataFrame) return crs end -function GeoInterface.geometrycolumns(table::DataFrame) +function GI.geometrycolumns(table::DataFrame) gc = DataAPI.metadata(table, "GEOINTERFACE:geometrycolumns", nothing) if isnothing(gc) # fall back to searching for "geometrycolumns" as a string gc = DataAPI.metadata(table, "geometrycolumns", (:geometry,)) @@ -75,29 +75,32 @@ end # feature interface, for use in generic code. And dispatch can always # handle a DataFrame by fixing the trait in a specialized method. +# GI.isfeaturecollection(::Type{<:DataFrame}) = true + # Here, we define a feature as a DataFrame row. -function GeoInterface.getfeature(df::DataFrame, i::Integer) +function GI.getfeature(df::DataFrame, i::Integer) return view(df, i, :) end # This is simply an optimized method, since we know what we have to do already. -GeoInterface.getfeature(df::DataFrame) = eachrow(df) +GI.getfeature(df::DataFrame) = eachrow(df) +# GI.isfeature(::Type{<:DataFrameRow}) = true # The geometry is defined as the first of the geometry columns. # TODO: how should we choose between the geometry columns? -function GeoInterface.geometry(row::DataFrameRow) - return row[first(GeoInterface.geometrycolumns(row))] +function GI.geometry(row::DataFrameRow) + return row[first(GI.geometrycolumns(row))] end # The properties are all other columns. -function GeoInterface.properties(row::DataFrameRow) - return row[DataFrames.Not(first(GeoInterface.geometrycolumns(row)))] +function GI.properties(row::DataFrameRow) + return row[DataFrames.Not(first(GI.geometrycolumns(row)))] end # Since `DataFrameRow` is simply a view of a DataFrame, we can reach back # to the original DataFrame to get the metadata. -GeoInterface.geometrycolumns(row::DataFrameRow) = - GeoInterface.geometrycolumns(getfield(row, :df)) # get the parent of the row view -GeoInterface.crs(row::DataFrameRow) = GeoInterface.crs(getfield(row, :df)) # get the parent of the row view +GI.geometrycolumns(row::DataFrameRow) = + GI.geometrycolumns(getfield(row, :df)) # get the parent of the row view +GI.crs(row::DataFrameRow) = GI.crs(getfield(row, :df)) # get the parent of the row view """ reproject(df::DataFrame, to_crs) @@ -146,7 +149,7 @@ end function reproject(sv::AbstractVector{<:AG.IGeometry}, from_crs, to_crs) Base.depwarn( - "`reproject(sv::AbstractVector)` will be deprecated in a future release." * + "`reproject(sv::AbstractVector)` will be deprecated in a future release. " * "Please use `reproject(df::DataFrame)` instead to make sure the dataframe crs metadata is updated.", :reproject, ) diff --git a/test/runtests.jl b/test/runtests.jl index 9eaa4d3..625969c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -15,33 +15,44 @@ REPO_URL = "https://github.com/yeesian/ArchGDALDatasets/blob/master/" remotefiles = [ ( - "ospy/data1/sites.dbf", + "https://github.com/yeesian/ArchGDALDatasets/blob/master/ospy/data1/sites.dbf?raw=true", "7df95edea06c46418287ae3430887f44f9116b29715783f7d1a11b2b931d6e7d", ), ( - "ospy/data1/sites.prj", + "https://github.com/yeesian/ArchGDALDatasets/blob/master/ospy/data1/sites.prj?raw=true", "81fb1a246728609a446b25b0df9ede41c3e7b6a133ce78f10edbd2647fc38ce1", ), ( - "ospy/data1/sites.sbn", + "https://github.com/yeesian/ArchGDALDatasets/blob/master/ospy/data1/sites.sbn?raw=true", "198d9d695f3e7a0a0ac0ebfd6afbe044b78db3e685fffd241a32396e8b341ed3", ), ( - "ospy/data1/sites.sbx", + "https://github.com/yeesian/ArchGDALDatasets/blob/master/ospy/data1/sites.sbx?raw=true", "49bbe1942b899d52cf1d1b01ea10bd481ec40bdc4c94ff866aece5e81f2261f6", ), ( - "ospy/data1/sites.shp", + "https://github.com/yeesian/ArchGDALDatasets/blob/master/ospy/data1/sites.shp?raw=true", "69af5a6184053f0b71f266dc54c944f1ec02013fb66dbb33412d8b1976d5ea2b", ), ( - "ospy/data1/sites.shx", + "https://github.com/yeesian/ArchGDALDatasets/blob/master/ospy/data1/sites.shx?raw=true", "1f3da459ccb151958743171e41e6a01810b2a007305d55666e01d680da7bbf08", ), + ( + "https://github.com/opengeospatial/geoparquet/raw/v1.0.0/examples/example.parquet", + "3dc1a3df76290cc62aa6d7aa6aa00d0988b3157ee77a042167b3f8302d05aa6c", + ), + ( + "https://storage.googleapis.com/open-geodata/linz-examples/nz-buildings-outlines.parquet", + "b64389237b9879c275aec4e47f0e6be8fdd5d64437a05aef10839f574efc5dbc", + ), + ( + "https://github.com/bjornharrtell/flatgeobuf/blob/master/test/data/countries.fgb?raw=true", + "d8dc3baf855320d10c6a662bf1171273849dd8a0935066b5b7b8dd83b3484cb3", + ), ] -for (f, sha) in remotefiles - localfn = joinpath(testdatadir, basename(f)) - url = REPO_URL * f * "?raw=true" +for (url, sha) in remotefiles + localfn = joinpath(testdatadir, basename(replace(url, "?raw=true" => ""))) PlatformEngines.download_verify(url, sha, localfn; force = true, quiet_download = false) end @@ -289,6 +300,64 @@ unknown_crs = GFT.WellKnownText( @test_throws ErrorException("Unable to open $wfn.") GDF.read(wfn) end + @testset "Shapefile" begin + using Shapefile + fn = joinpath(testdatadir, "sites.shp") + df = GDF.read(fn) + df2 = GDF.read(GDF.ArchGDALDriver(), fn) + @test names(df) == names(df2) + @test nrow(df) == nrow(df2) + @test GI.trait(df.geometry[1]) == GI.trait(df2.geometry[1]) + @test GI.coordinates(df.geometry[1]) == GI.coordinates(df2.geometry[1]) + + GDF.write("test_native.shp", df; force = true) + GDF.write(GDF.ArchGDALDriver(), "test.shp", df; force = true) + end + @testset "GeoJSON" begin + using GeoJSON + fn = joinpath(testdatadir, "test_polygons.geojson") + df = GDF.read(fn) + df2 = GDF.read(GDF.ArchGDALDriver(), fn) + @test names(df) == names(df2) + @test nrow(df) == nrow(df2) + @test GI.trait(df.geometry[1]) == GI.trait(df2.geometry[1]) + @test all( + isapprox.( + collect.(GI.coordinates(df.geometry[1])[1]), + GI.coordinates(df2.geometry[1])[1], + ), + ) + GDF.write("test_native.geojson", df) + GDF.write(GDF.ArchGDALDriver(), "test.geojson", df) + end + @testset "FlatGeobuf" begin + using FlatGeobuf + fn = joinpath(testdatadir, "countries.fgb") + df = GDF.read(fn) + df2 = GDF.read(GDF.ArchGDALDriver(), fn) + @test sort(names(df)) == sort(names(df2)) + @test nrow(df) == nrow(df2) + # FlatGeobuf does not support GeoInterface yet + @test_broken GI.trait(df.geometry[1]) == GI.trait(df2.geometry[1]) + @test_broken GI.coordinates(df.geometry[1]) == GI.coordinates(df2.geometry[1]) + + # GDF.write("test_native.fgb", df) # Can't convert FlatGeobuf to ArchGDAL + GDF.write("test_native.fgb", df2) + GDF.write(GDF.ArchGDALDriver(), "test.fgb", df2) + end + @testset "GeoParquet" begin + using GeoParquet + fn = joinpath(testdatadir, "example.parquet") + df = GDF.read(fn) + df2 = GDF.read(GDF.ArchGDALDriver(), fn) + @test sort(names(df)) == sort(names(df2)) + @test nrow(df) == nrow(df2) + @test GI.trait(df.geometry[1]) == GI.trait(df2.geometry[1]) + @test GI.coordinates(df.geometry[1]) == GI.coordinates(df2.geometry[1]) + + GDF.write("test_native.parquet", df) + GDF.write(GDF.ArchGDALDriver(), "test.parquet", df) + end @testset "Writing crs of geometry" begin geom = GI.Wrappers.Point(0, 0; crs = GFT.EPSG(4326)) df = DataFrame(; geometry = [geom])