From 2a8c80c125b2f46431bf3875a34287bb6343d7c4 Mon Sep 17 00:00:00 2001 From: Maarten Pronk Date: Sat, 23 Nov 2024 17:33:10 +0100 Subject: [PATCH] Add on table, and derive crs from geometry. --- README.md | 4 +++ docs/src/tutorials/ops.md | 6 +++- src/GeoDataFrames.jl | 3 -- src/exports.jl | 38 ++++++++++++++++++++- src/utils.jl | 72 ++++++++++++++++++++++++++++++++++++++- test/runtests.jl | 20 +++++++++-- 6 files changed, 134 insertions(+), 9 deletions(-) diff --git a/README.md b/README.md index 9096ae7..896255d 100644 --- a/README.md +++ b/README.md @@ -96,6 +96,10 @@ df ### Reprojection ```julia +#Either reproject the whole dataframe (which will set the correct metadata): +dfr = reproject(df, EPSG(4326), EPSG(28992)) + +# Or reproject the individual geometries only: df.geometry = reproject(df.geometry, EPSG(4326), EPSG(28992)) 10-element Vector{ArchGDAL.IGeometry{ArchGDAL.wkbPolygon}}: Geometry: POLYGON ((-472026.042542408 -4406233.59953401,-537 ... 401)) diff --git a/docs/src/tutorials/ops.md b/docs/src/tutorials/ops.md index 9b75be0..5b23a73 100644 --- a/docs/src/tutorials/ops.md +++ b/docs/src/tutorials/ops.md @@ -25,7 +25,11 @@ df ## Reprojection ```julia -df.geom = reproject(df.geom, EPSG(4326), EPSG(28992)) +#Either reproject the whole dataframe (which will set the correct metadata): +dfr = reproject(df, EPSG(4326), EPSG(28992)) + +# Or reproject the individual geometries only: +df.geometry = reproject(df.geometry, EPSG(4326), EPSG(28992)) 10-element Vector{ArchGDAL.IGeometry{ArchGDAL.wkbPolygon}}: Geometry: POLYGON ((-472026.042542408 -4406233.59953401,-537 ... 401)) Geometry: POLYGON ((-417143.506054105 -4395423.99277048,-482 ... 048)) diff --git a/src/GeoDataFrames.jl b/src/GeoDataFrames.jl index 40ea4d1..b5fccdd 100644 --- a/src/GeoDataFrames.jl +++ b/src/GeoDataFrames.jl @@ -8,9 +8,6 @@ import GeoInterface using DataAPI using Reexport -@reexport using Extents -@reexport using GeoFormatTypes - include("exports.jl") include("io.jl") include("utils.jl") diff --git a/src/exports.jl b/src/exports.jl index f379b2d..ea040db 100644 --- a/src/exports.jl +++ b/src/exports.jl @@ -11,7 +11,43 @@ createpolygon, createmultilinestring, createmultipolygon -@reexport using ArchGDAL: reproject + +@reexport using Extents +@reexport using GeoInterface: crs + +using GeoFormatTypes: + AbstractWellKnownText, + CoordSys, + CoordinateReferenceSystemFormat, + EPSG, + ESRIWellKnownText, + GML, + GeoFormat, + GeoFormatTypes, + GeometryFormat, + KML, + MixedFormat, + ProjJSON, + ProjString, + WellKnownBinary, + WellKnownText, + WellKnownText2 #=GeoJSON,=# +export AbstractWellKnownText, + CoordSys, + CoordinateReferenceSystemFormat, + EPSG, + ESRIWellKnownText, + GML, + GeoFormat, + GeoFormatTypes, + GeometryFormat, + KML, + MixedFormat, + ProjJSON, + ProjString, + WellKnownBinary, + WellKnownText, + WellKnownText2 #=GeoJSON,=# AG.intersects(a::Vector{AG.IGeometry{T}}, b::Vector{AG.IGeometry{X}}) where {X, T} = AG.intersects.(a, b) diff --git a/src/utils.jl b/src/utils.jl index df82781..7ad7e5d 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -23,12 +23,26 @@ end function getcrs(table) if GeoInterface.isfeaturecollection(table) return GeoInterface.crs(table) - elseif first(DataAPI.metadatasupport(typeof(table))) + end + crs = geomcrs(table) + if !isnothing(crs) + return crs + end + if first(DataAPI.metadatasupport(typeof(table))) 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) end return crs + end + nothing +end + +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)) else return nothing end @@ -84,3 +98,59 @@ end 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 + +""" + reproject(df::DataFrame, to_crs) + +Reproject the geometries in a DataFrame `df` to a new Coordinate Reference System `to_crs`, from the current CRS. +See also [`reproject(df, from_crs, to_crs)`](@ref) and the in place version [`reproject!(df, to_crs)`](@ref). +""" +function reproject(df::DataFrame, to_crs) + reproject!(copy(df), to_crs) +end + +""" + reproject(df::DataFrame, from_crs, to_crs) + +Reproject the geometries in a DataFrame `df` from the crs `from_crs` to a new crs `to_crs`. +This overrides any current CRS of the Dataframe. +""" +function reproject(df::DataFrame, from_crs, to_crs) + reproject!(copy(df), from_crs, to_crs) +end + +""" + reproject!(df::DataFrame, to_crs) + +Reproject the geometries in a DataFrame `df` to a new Coordinate Reference System `to_crs`, from the current CRS, in place. +""" +function reproject!(df::DataFrame, to_crs) + reproject!(df, getcrs(df), to_crs) +end + +""" + reproject!(df::DataFrame, from_crs, to_crs) + +Reproject the geometries in a DataFrame `df` from the crs `from_crs` to a new crs `to_crs` in place. +This overrides any current CRS of the Dataframe. +""" +function reproject!(df::DataFrame, from_crs, to_crs) + columns = Tables.columns(df) + for gc in getgeometrycolumns(df) + gc in Tables.columnnames(columns) || continue + reproject(df[!, gc], from_crs, to_crs) + end + metadata!(df, "crs", to_crs; style = :note) + metadata!(df, "GEOINTERFACE:crs", to_crs; style = :note) +end + +function reproject(sv::AbstractVector{<:AG.IGeometry}, from_crs, to_crs) + Base.depwarn( + "`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, + ) + AG.reproject.(sv, Ref(from_crs), Ref(to_crs)) +end + +export reproject, reproject! diff --git a/test/runtests.jl b/test/runtests.jl index 9fc3ad8..9eaa4d3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -192,8 +192,10 @@ unknown_crs = GFT.WellKnownText( @testset "Reproject" begin 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 + geoms = GDF.reproject(AG.clone.(table.geometry), GFT.EPSG(4326), GFT.EPSG(28992)) + ntable = GDF.reproject(table, GFT.EPSG(4326), GFT.EPSG(28992)) + @test GDF.AG.getpoint(geoms[1], 0)[1] ≈ -587791.596556932 + @test GDF.AG.getpoint(ntable.geometry[1], 0)[1] ≈ -587791.596556932 GDF.write( joinpath(testdatadir, "test_reprojection.gpkg"), table; @@ -269,11 +271,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") + isdir(gdbdir) && 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 +288,16 @@ 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 "Writing crs of geometry" begin + geom = GI.Wrappers.Point(0, 0; crs = GFT.EPSG(4326)) + df = DataFrame(; geometry = [geom]) + @test isnothing(GI.crs(df)) + GDF.write("test_geom_crs.gpkg", df) + df = GDF.read("test_geom_crs.gpkg") + @test GI.crs(df) == GFT.WellKnownText{GFT.CRS}( + GFT.CRS(), + "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]", + ) + end end