From 466b7a8906c9bbea6837520f9eb35b15b57c87da Mon Sep 17 00:00:00 2001 From: Maarten Pronk Date: Wed, 20 Nov 2024 10:49:33 +0100 Subject: [PATCH 1/3] Add native driver extensions. --- Project.toml | 18 ++++++- src/GeoDataFrames.jl | 6 ++- src/io.jl | 119 ++++++++++++++++++++++++++----------------- src/utils.jl | 32 ++++++------ test/runtests.jl | 90 ++++++++++++++++++++++++++++---- 5 files changed, 191 insertions(+), 74 deletions(-) diff --git a/Project.toml b/Project.toml index df95f29..4e46de3 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] +GeoJSON = "61d90e0f-e114-555e-ac52-39dfb47a3ef9" +GeoParquet = "e99870d8-ce00-4fdd-aeee-e09192881159" +Shapefile = "8e980c4a-a4fe-5da2-b3a7-4b4b0353a2f4" +FlatGeobuf = "d985ece1-97de-4d33-914c-38fb84042e15" + +[extensions] +GeoJSONExt = "GeoJSON" +GeoParquetExt = "GeoParquet" +ShapefileExt = "Shapefile" +FlatGeobufExt = "FlatGeobuf" + [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.6" @@ -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/src/GeoDataFrames.jl b/src/GeoDataFrames.jl index 40ea4d1..600fdf4 100644 --- a/src/GeoDataFrames.jl +++ b/src/GeoDataFrames.jl @@ -4,7 +4,7 @@ import ArchGDAL as AG using DataFrames using Tables import GeoFormatTypes as GFT -import GeoInterface +import GeoInterface as GI using DataAPI using Reexport @@ -12,7 +12,11 @@ using Reexport @reexport using GeoFormatTypes include("exports.jl") +include("sources.jl") include("io.jl") include("utils.jl") +function load end +function save end + end # module diff --git a/src/io.jl b/src/io.jl index 1aa736c..97e211d 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,48 @@ 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; kwargs...) startswith(fn, "/vsi") || occursin(":", fn) || isfile(fn) || isdir(fn) || error("File not found.") + + if :layer in keys(kwargs) + layer = kwargs[:layer] + kwargs = (k => v for (k, v) in pairs(kwargs) if k != :layer) + else + layer = nothing + end + 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.") - 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( @@ -110,7 +123,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 +150,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 +158,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 +186,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 +284,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 df82781..1154a67 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) elseif first(DataAPI.metadatasupport(typeof(table))) crs = DataAPI.metadata(table, "GEOINTERFACE:crs", nothing) if isnothing(crs) # fall back to searching for "crs" as a string @@ -38,7 +38,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) @@ -46,7 +46,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,)) @@ -61,26 +61,28 @@ 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 diff --git a/test/runtests.jl b/test/runtests.jl index 9fc3ad8..9c51484 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 @@ -269,11 +280,11 @@ unknown_crs = GFT.WellKnownText( @testset "Read geodatabase (folder)" begin table = DataFrame(; geom = AG.createpoint(1.0, 2.0), name = "test") gdbdir = joinpath(testdatadir, "test_options.gdb") + rm(gdbdir; recursive = true) GDF.write(gdbdir, table; driver = "OpenFileGDB", geom_column = :geom) @test isdir(gdbdir) table = GDF.read(gdbdir) @test nrow(table) == 1 - rm(gdbdir; recursive = true) end @testset "Non-spatial columns #77" begin @@ -286,4 +297,63 @@ unknown_crs = GFT.WellKnownText( wfn = "C:\\non_existing_folder\\non_existing_file.shp" @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 end From c522d269622c8ebc796f70364314ec8aa6169389 Mon Sep 17 00:00:00 2001 From: Maarten Pronk Date: Wed, 20 Nov 2024 11:51:44 +0100 Subject: [PATCH 2/3] Let's also add files. --- ext/FlatGeobufExt.jl | 15 +++++++++++++++ ext/GeoJSONExt.jl | 13 +++++++++++++ ext/GeoParquetExt.jl | 14 ++++++++++++++ ext/ShapefileExt.jl | 13 +++++++++++++ src/sources.jl | 33 +++++++++++++++++++++++++++++++++ 5 files changed, 88 insertions(+) create mode 100644 ext/FlatGeobufExt.jl create mode 100644 ext/GeoJSONExt.jl create mode 100644 ext/GeoParquetExt.jl create mode 100644 ext/ShapefileExt.jl create mode 100644 src/sources.jl diff --git a/ext/FlatGeobufExt.jl b/ext/FlatGeobufExt.jl new file mode 100644 index 0000000..abaf14a --- /dev/null +++ b/ext/FlatGeobufExt.jl @@ -0,0 +1,15 @@ +module FlatGeobufExt +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/GeoJSONExt.jl b/ext/GeoJSONExt.jl new file mode 100644 index 0000000..932a2d6 --- /dev/null +++ b/ext/GeoJSONExt.jl @@ -0,0 +1,13 @@ +module GeoJSONExt +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/GeoParquetExt.jl b/ext/GeoParquetExt.jl new file mode 100644 index 0000000..f72f61a --- /dev/null +++ b/ext/GeoParquetExt.jl @@ -0,0 +1,14 @@ +module GeoParquetExt + +using GeoDataFrames: GeoParquetDriver, GeoDataFrames +import GeoInterface as GI +using GeoParquet + +function GeoDataFrames.read(::GeoParquetDriver, fname::AbstractString; kwargs...) + GeoParquet.read(fname; kwargs...) +end + +function GeoDataFrames.write(::GeoParquetDriver, fname::AbstractString, data; kwargs...) + GeoParquet.write(fname, data, GI.geometrycolumns(data); kwargs...) +end +end diff --git a/ext/ShapefileExt.jl b/ext/ShapefileExt.jl new file mode 100644 index 0000000..ae8869b --- /dev/null +++ b/ext/ShapefileExt.jl @@ -0,0 +1,13 @@ +module ShapefileExt +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/sources.jl b/src/sources.jl new file mode 100644 index 0000000..53b3bf6 --- /dev/null +++ b/src/sources.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" From 07b144e6fcf51388e070ce3b96b1cb1e91336dbf Mon Sep 17 00:00:00 2001 From: Maarten Pronk Date: Tue, 7 Jan 2025 23:27:18 -0100 Subject: [PATCH 3/3] Rename geom to geometry. Address comments. --- src/GeoDataFrames.jl | 2 +- src/{sources.jl => drivers.jl} | 0 src/io.jl | 13 +++++-------- 3 files changed, 6 insertions(+), 9 deletions(-) rename src/{sources.jl => drivers.jl} (100%) diff --git a/src/GeoDataFrames.jl b/src/GeoDataFrames.jl index 600fdf4..d857351 100644 --- a/src/GeoDataFrames.jl +++ b/src/GeoDataFrames.jl @@ -12,7 +12,7 @@ using Reexport @reexport using GeoFormatTypes include("exports.jl") -include("sources.jl") +include("drivers.jl") include("io.jl") include("utils.jl") diff --git a/src/sources.jl b/src/drivers.jl similarity index 100% rename from src/sources.jl rename to src/drivers.jl diff --git a/src/io.jl b/src/io.jl index 97e211d..f39f68a 100644 --- a/src/io.jl +++ b/src/io.jl @@ -60,20 +60,13 @@ function read(driver::AbstractDriver, fn::AbstractString; kwargs...) read(ArchGDALDriver(), fn; kwargs...) end -function read(driver::ArchGDALDriver, fn::AbstractString; kwargs...) +function read(driver::ArchGDALDriver, fn::AbstractString; layer=nothing, kwargs...) startswith(fn, "/vsi") || occursin(":", fn) || isfile(fn) || isdir(fn) || error("File not found.") - if :layer in keys(kwargs) - layer = kwargs[:layer] - kwargs = (k => v for (k, v) in pairs(kwargs) if k != :layer) - else - layer = nothing - end - t = AG.read(fn; kwargs...) do ds ds.ptr == C_NULL && error("Unable to open $fn.") if AG.nlayer(ds) > 1 && isnothing(layer) @@ -107,6 +100,10 @@ function read(::ArchGDALDriver, ds, layer) rename!(df, Symbol("") => :geometry) replace!(gnames, Symbol("") => :geometry) end + if "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) metadata!(df, "crs", crs; style = :note)