From d183b95d90073274c2be5b6191b5688c976d32fe Mon Sep 17 00:00:00 2001
From: Elias Carvalho <eliascarvdev@gmail.com>
Date: Tue, 1 Oct 2024 13:50:16 -0300
Subject: [PATCH 1/4] Add support for saving CRS in Shapefile

---
 Project.toml         |  2 +-
 src/conversion.jl    | 11 +++++++++++
 src/utils.jl         | 11 +++++++++++
 test/gis.jl          |  6 +++---
 test/io/shapefile.jl |  8 ++++----
 5 files changed, 30 insertions(+), 8 deletions(-)

diff --git a/Project.toml b/Project.toml
index 4488b69..5f61674 100644
--- a/Project.toml
+++ b/Project.toml
@@ -39,7 +39,7 @@ ArchGDAL = "0.10"
 CSV = "0.10"
 Colors = "0.12"
 CommonDataModel = "0.2, 0.3"
-CoordRefSystems = "0.12 - 0.14"
+CoordRefSystems = "0.15"
 FileIO = "1.16"
 Format = "1.3"
 GRIBDatasets = "0.3"
diff --git a/src/conversion.jl b/src/conversion.jl
index 11ac5a7..0a7d223 100644
--- a/src/conversion.jl
+++ b/src/conversion.jl
@@ -17,6 +17,8 @@ GI.geomtrait(::MultiRope) = GI.MultiLineStringTrait()
 GI.geomtrait(::MultiRing) = GI.MultiLineStringTrait()
 GI.geomtrait(::MultiPolygon) = GI.MultiPolygonTrait()
 
+GI.crs(geom::Geometry) = gicrs(crs(geom))
+
 GI.ncoord(::GI.PointTrait, p::Point) = CoordRefSystems.ncoords(crs(p))
 GI.getcoord(::GI.PointTrait, p::Point) = ustrip.(raw(coords(p)))
 GI.getcoord(trait::GI.PointTrait, p::Point, i) = GI.getcoord(trait, p)[i]
@@ -42,9 +44,18 @@ GI.getgeom(::GI.AbstractGeometryTrait, m::Multi, i) = parent(m)[i]
 
 GI.isfeaturecollection(::Type{<:AbstractGeoTable}) = true
 GI.trait(::AbstractGeoTable) = GI.FeatureCollectionTrait()
+GI.crs(gtb::AbstractGeoTable) = gicrs(crs(domain(gtb)))
 GI.nfeature(::Any, gtb::AbstractGeoTable) = nrow(gtb)
 GI.getfeature(::Any, gtb::AbstractGeoTable, i) = gtb[i, :]
 
+function gicrs(CRS)
+  try
+    GFT.WellKnownText2(GFT.CRS(), wktstring(CoordRefSystems.code(CRS)))
+  catch
+    nothing
+  end
+end
+
 # --------------------------------------
 # Convert geometries to Meshes.jl types
 # --------------------------------------
diff --git a/src/utils.jl b/src/utils.jl
index f851987..c6291d3 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -47,3 +47,14 @@ function uniquenames(names, newnames)
     uniquename(names, name)
   end
 end
+
+function wktstring(code; format="WKT2", multiline=false)
+  spref = AG.importUserInput(codestring(code))
+  options = ["FORMAT=$format", "MULTILINE=$(multiline ? "YES" : "NO")"]
+  wktptr = Ref{Cstring}()
+  AG.GDAL.osrexporttowktex(spref, wktptr, options)
+  unsafe_string(wktptr[])
+end
+
+codestring(::Type{EPSG{Code}}) where {Code} = "EPSG:$Code"
+codestring(::Type{ESRI{Code}}) where {Code} = "ESRI:$Code"
diff --git a/test/gis.jl b/test/gis.jl
index 65309fb..c5c317c 100644
--- a/test/gis.jl
+++ b/test/gis.jl
@@ -20,21 +20,21 @@
   file = joinpath(savedir, "gis-points.shp")
   GeoIO.save(file, gtpoint)
   gtb = GeoIO.load(file)
-  @test_broken gtb.geometry == gtpoint.geometry
+  @test gtb.geometry == gtpoint.geometry
   @test values(gtb) == values(gtpoint)
 
   # note: Shapefile saves Chain as MultiChain
   file = joinpath(savedir, "gis-rings.shp")
   GeoIO.save(file, gtring)
   gtb = GeoIO.load(file)
-  @test_broken _isequal(gtb.geometry, gtring.geometry)
+  @test _isequal(gtb.geometry, gtring.geometry)
   @test values(gtb) == values(gtring)
 
   # note: Shapefile saves PolyArea as MultiPolyArea
   file = joinpath(savedir, "gis-polys.shp")
   GeoIO.save(file, gtpoly)
   gtb = GeoIO.load(file)
-  @test_broken _isequal(gtb.geometry, gtpoly.geometry)
+  @test _isequal(gtb.geometry, gtpoly.geometry)
   @test values(gtb) == values(gtpoly)
 
   # GeoJSON
diff --git a/test/io/shapefile.jl b/test/io/shapefile.jl
index 9abdf76..3cb6b27 100644
--- a/test/io/shapefile.jl
+++ b/test/io/shapefile.jl
@@ -68,9 +68,9 @@
     file1 = joinpath(datadir, "points.shp")
     file2 = joinpath(savedir, "points.shp")
     gtb1 = GeoIO.load(file1)
-    GeoIO.save(file2, gtb1)
+    GeoIO.save(file2, gtb1, force=true)
     gtb2 = GeoIO.load(file2)
-    @test_broken gtb1.geometry == gtb2.geometry
+    @test gtb1.geometry == gtb2.geometry
     @test values(gtb1) == values(gtb2)
     @test values(gtb1, 0) == values(gtb2, 0)
 
@@ -79,7 +79,7 @@
     gtb1 = GeoIO.load(file1)
     GeoIO.save(file2, gtb1)
     gtb2 = GeoIO.load(file2)
-    @test_broken gtb1.geometry == gtb2.geometry
+    @test gtb1.geometry == gtb2.geometry
     @test values(gtb1) == values(gtb2)
     @test values(gtb1, 0) == values(gtb2, 0)
 
@@ -88,7 +88,7 @@
     gtb1 = GeoIO.load(file1)
     GeoIO.save(file2, gtb1)
     gtb2 = GeoIO.load(file2)
-    @test_broken gtb1.geometry == gtb2.geometry
+    @test gtb1.geometry == gtb2.geometry
     @test values(gtb1) == values(gtb2)
     @test values(gtb1, 0) == values(gtb2, 0)
   end

From 0ca085c35ebaa13b15cc7841fa7a30146a6a1579 Mon Sep 17 00:00:00 2001
From: Elias Carvalho <eliascarvdev@gmail.com>
Date: Tue, 1 Oct 2024 15:13:00 -0300
Subject: [PATCH 2/4] Fix typo

---
 test/io/shapefile.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/io/shapefile.jl b/test/io/shapefile.jl
index 3cb6b27..7b4073d 100644
--- a/test/io/shapefile.jl
+++ b/test/io/shapefile.jl
@@ -68,7 +68,7 @@
     file1 = joinpath(datadir, "points.shp")
     file2 = joinpath(savedir, "points.shp")
     gtb1 = GeoIO.load(file1)
-    GeoIO.save(file2, gtb1, force=true)
+    GeoIO.save(file2, gtb1)
     gtb2 = GeoIO.load(file2)
     @test gtb1.geometry == gtb2.geometry
     @test values(gtb1) == values(gtb2)

From 400d935c309f7e58270e9a59b6369bc7aa5082db Mon Sep 17 00:00:00 2001
From: Elias Carvalho <eliascarvdev@gmail.com>
Date: Tue, 1 Oct 2024 17:04:40 -0300
Subject: [PATCH 3/4] Update helper functions

---
 src/utils.jl | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/utils.jl b/src/utils.jl
index c6291d3..2e46d15 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -49,12 +49,14 @@ function uniquenames(names, newnames)
 end
 
 function wktstring(code; format="WKT2", multiline=false)
-  spref = AG.importUserInput(codestring(code))
-  options = ["FORMAT=$format", "MULTILINE=$(multiline ? "YES" : "NO")"]
+  spref = spatialref(code)
   wktptr = Ref{Cstring}()
+  options = ["FORMAT=$format", "MULTILINE=$(multiline ? "YES" : "NO")"]
   AG.GDAL.osrexporttowktex(spref, wktptr, options)
   unsafe_string(wktptr[])
 end
 
+spatialref(code) = AG.importUserInput(codestring(code))
+
 codestring(::Type{EPSG{Code}}) where {Code} = "EPSG:$Code"
 codestring(::Type{ESRI{Code}}) where {Code} = "ESRI:$Code"

From a759d93eca16b030e4348e3d4d0ad71eb356a12f Mon Sep 17 00:00:00 2001
From: Elias Carvalho <eliascarvdev@gmail.com>
Date: Tue, 1 Oct 2024 17:21:08 -0300
Subject: [PATCH 4/4] Update code

---
 src/GeoIO.jl | 1 +
 src/utils.jl | 2 +-
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/GeoIO.jl b/src/GeoIO.jl
index e7a46f9..6e46a69 100644
--- a/src/GeoIO.jl
+++ b/src/GeoIO.jl
@@ -49,6 +49,7 @@ import ArchGDAL as AG
 import GeoParquet as GPQ
 import GeoInterface as GI
 import GeoFormatTypes as GFT
+import ArchGDAL.GDAL
 
 # image extensions
 const IMGEXTS = [".png", ".jpg", ".jpeg"]
diff --git a/src/utils.jl b/src/utils.jl
index 2e46d15..060ab0e 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -52,7 +52,7 @@ function wktstring(code; format="WKT2", multiline=false)
   spref = spatialref(code)
   wktptr = Ref{Cstring}()
   options = ["FORMAT=$format", "MULTILINE=$(multiline ? "YES" : "NO")"]
-  AG.GDAL.osrexporttowktex(spref, wktptr, options)
+  GDAL.osrexporttowktex(spref, wktptr, options)
   unsafe_string(wktptr[])
 end