Add support for saving CRS in Shapefile

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"

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

src/conversion.jl

@@ -17,6 +17,8 @@
 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.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
 # --------------------------------------

src/utils.jl

@@ -47,3 +47,16 @@
     uniquename(names, name)
   end
 end
+
+function wktstring(code; format="WKT2", multiline=false)
+  spref = spatialref(code)
+  wktptr = Ref{Cstring}()
+  options = ["FORMAT=$format", "MULTILINE=$(multiline ? "YES" : "NO")"]
+  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"

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

test/io/shapefile.jl

@@ -70,7 +70,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)
 
@@ -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