Add support for converting to WKT to EPSG (#391)

yeesian · Nov 28, 2023 · b6be106 · b6be106
1 parent 1671e1f
commit b6be106
Show file tree

Hide file tree

Showing 3 changed files with 38 additions and 0 deletions.
diff --git 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
@@ -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
@@ -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