Add reproject on table, and derive crs from geometry

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))

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))

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")

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)

src/utils.jl

@@ -23,12 +23,26 @@
 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 @@
 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!

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