From b6be1068655f2921159c32d7b8bb9c4d836f885f Mon Sep 17 00:00:00 2001
From: Alex Gardner <32276930+alex-s-gardner@users.noreply.github.com>
Date: Tue, 28 Nov 2023 00:41:45 -0800
Subject: [PATCH] Add support for converting to WKT to EPSG (#391)

---
 src/convert.jl       |  3 +++
 src/spatialref.jl    | 27 +++++++++++++++++++++++++++
 test/test_convert.jl |  8 ++++++++
 3 files changed, 38 insertions(+)

diff --git a/src/convert.jl b/src/convert.jl
index cff86554..9b1ddc27 100644
--- a/src/convert.jl
+++ b/src/convert.jl
@@ -112,6 +112,9 @@ end
 function unsafe_convertcrs(::Type{<:GFT.ProjString}, crsref)
     return GFT.ProjString(toPROJ4(crsref))
 end
+function unsafe_convertcrs(::Type{<:GFT.EPSG}, crsref)
+    return GFT.EPSG(toEPSG(crsref))
+end
 function unsafe_convertcrs(::Type{<:GFT.WellKnownText}, crsref)
     return GFT.WellKnownText(GFT.CRS(), toWKT(crsref))
 end
diff --git a/src/spatialref.jl b/src/spatialref.jl
index 0d1fef70..01957258 100644
--- a/src/spatialref.jl
+++ b/src/spatialref.jl
@@ -605,6 +605,7 @@ Convert this SRS into a nicely formatted WKT string for display to a person.
 function toWKT(spref::AbstractSpatialRef, simplify::Bool)::String
     wktptr = Ref{Cstring}()
     result = GDAL.osrexporttoprettywkt(spref, wktptr, simplify)
+
     @ogrerr result "Failed to convert this SRS into pretty WKT"
     return unsafe_string(wktptr[])
 end
@@ -621,6 +622,32 @@ function toPROJ4(spref::AbstractSpatialRef)::String
     return unsafe_string(projptr[])
 end
 
+"""
+    toEPSG(spref::AbstractSpatialRef)
+
+Export EPSG code for this coordinate system if available.
+"""
+function toEPSG(spref::AbstractSpatialRef)::Int64
+    result = GDAL.osrgetauthorityname(spref.ptr, "PROJCS")
+
+    if !isnothing(result)
+        projcs = "PROJCS"
+    else
+        result = GDAL.osrgetauthorityname(spref.ptr, "GEOGCS")
+        projcs = "GEOGCS"
+    end
+
+    if isnothing(result)
+        error("No PROJCS or GEOGCS Authority found")
+    elseif result == "EPSG"
+            epsg = GDAL.osrgetauthoritycode(spref.ptr, projcs)
+            epsg = parse(Int64, epsg)
+        return epsg
+    else
+        error("$result is not an EPSG authority")
+    end
+end
+
 """
     toXML(spref::AbstractSpatialRef)
 
diff --git a/test/test_convert.jl b/test/test_convert.jl
index 3310e9f0..327a8f2a 100644
--- a/test/test_convert.jl
+++ b/test/test_convert.jl
@@ -40,6 +40,14 @@ import GeoFormatTypes as GFT
         ) == proj4326
         @test convert(GFT.CoordSys, GFT.CRS(), proj4326) isa GFT.CoordSys
         @test convert(GFT.GML, GFT.CRS(), proj4326) isa GFT.GML
+
+        epsg = GFT.EPSG(4326)
+        wkt = convert(GFT.WellKnownText, epsg)
+        @test convert(GFT.EPSG, wkt) == epsg
+
+        epsg = GFT.EPSG(3013)
+        wkt = convert(GFT.WellKnownText, epsg)
+        @test convert(GFT.EPSG, wkt) == epsg
     end
 
     @testset "geometry conversions" begin