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 087f26e4..3af3f968 100644 --- a/src/spatialref.jl +++ b/src/spatialref.jl @@ -600,6 +600,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 @@ -616,6 +617,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