Add support for converting GeoReferenceTypes to EPSG

[alex-s-gardner]

currently ArchGDAL can not do the following operation: 'a = convert(GFT.EPSG, GFT.WellKnownText(spatial_ref))' [see end of post].

But it can convert between other GeoReferenceTypes refs e.g. a = convert(GFT.ProjString, GFT.WellKnownText(spatial_ref))

EPSG is a bit funny because there is no GDAL.osrexporttoepsg and instead the Authority needs to be parsed from WellKnownText. To support the conversion I think we just need to add 2 function: [1] unsafe_convertcrs(::Type{<:GFT.EPSG}, crsref) [2] toEPSG(spref::AbstractSpatialRef)::String

In a branch I added 'function unsafe_convertcrs(::Type{<:GFT.EPSG}, crsref)' to convert.jl without issue but when it came to adding function toEPSG(spref::AbstractSpatialRef)::String to spatialref.jl I couldn't figure out a robust way to extract the EPSG code from the WKT.

Does anyone have a suggestion for how toEPSG should be written?

wkt = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"
"GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"

julia> a = convert(GFT.EPSG, GFT.WellKnownText(wkt))
ERROR: MethodError: no method matching unsafe_convertcrs(::Type{GeoFormatTypes.EPSG}, ::ArchGDAL.ISpatialRef)

Closest candidates are:
  unsafe_convertcrs(::Type{<:GeoFormatTypes.CoordSys}, ::Any)
   @ ArchGDAL ~/.julia/packages/ArchGDAL/AZONk/src/convert.jl:109
  unsafe_convertcrs(::Type{<:GeoFormatTypes.ProjString}, ::Any)
   @ ArchGDAL ~/.julia/packages/ArchGDAL/AZONk/src/convert.jl:112
  unsafe_convertcrs(::Type{<:GeoFormatTypes.WellKnownText}, ::Any)
   @ ArchGDAL ~/.julia/packages/ArchGDAL/AZONk/src/convert.jl:115
  ...

Stacktrace:
 [1] convert(target::Type{GeoFormatTypes.EPSG}, mode::GeoFormatTypes.CRS, source::GeoFormatTypes.WellKnownText{GeoFormatTypes.Unknown})
   @ ArchGDAL ~/.julia/packages/ArchGDAL/AZONk/src/convert.jl:106
 [2] convert(target::Type{GeoFormatTypes.EPSG}, source::GeoFormatTypes.WellKnownText{GeoFormatTypes.Unknown}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ GeoFormatTypes ~/.julia/packages/GeoFormatTypes/Lkriv/src/GeoFormatTypes.jl:118
 [3] convert(target::Type{GeoFormatTypes.EPSG}, source::GeoFormatTypes.WellKnownText{GeoFormatTypes.Unknown})
   @ GeoFormatTypes ~/.julia/packages/GeoFormatTypes/Lkriv/src/GeoFormatTypes.jl:95
 [4] top-level scope
   @ REPL[26]:1

[rafaqz]

julia> a = convert(GFT.EPSG, GFT.WellKnownText(spatial_ref))

I dont get it, EPSG doesnt support spatial refs... its just codes for projections right?

GeoFormatTypes/ArchGDAL only have the spatial ref conversions for types that can hold both, like WKT and WNB

You also define wkt in the example but dont use it?

[alex-s-gardner]

@rafaqz I want to be able to do the inverse of what's shown in the GeoFormatTypes.jl ReadMe example:

convert(WellKnownText, EPSG(4326))

does that make sense?

[rafaqz]

Oh right I forgot the terminology. But best if you use a full MWE so we can follow whats happening right through, I'm not sure what spatial_ref is in your code.

But no I don't know if GDAL can do that, seems I didn't find it when I implemented these things originally.

Doing it manually seems like it wont be robust. I'm not sure authority is required or that it has to be EPSG.

[visr]

In PROJ there is this: https://proj.org/en/9.2/development/reference/functions.html#c.proj_get_id_code

Don't know if something like that is also availabe in the GDAL C API though.

[alex-s-gardner]

I'm not sure what spatial_ref is in your code

That was a mistake on my part... should have been wkt. Now updated

[alex-s-gardner]

I found GDAL.osrgetauthoritycode

which I think should be used in function [2] from op as:

function toEPSG(spref::AbstractSpatialRef)::String
    epsg = Ref{Cstring}()
    result = GDAL.osrgetauthoritycode(spref, epsg)
    @ogrerr result "Failed to retrieve EPSG code for this SRS"
    return unsafe_string(epsg[])
end

When I implement [1] and [2] in ArchGDAL I now get cbinding errors that are over my head:

convert(GFT.EPSG, wkt)

Closest candidates are:
  convert(::Type{Cstring}, ::Union{Ptr{Int8}, Ptr{Nothing}, Ptr{UInt8}})
   @ Base c.jl:167
  convert(::Type{T}, ::T) where T
   @ Base Base.jl:64

Stacktrace:
 [1] cconvert(T::Type, x::Base.RefValue{Cstring})
   @ Base ./essentials.jl:492
 [2] osrgetauthoritycode(hSRS::ArchGDAL.ISpatialRef, pszTargetKey::Base.RefValue{Cstring})
   @ GDAL ~/.julia/packages/GDAL/A0qjS/src/libgdal.jl:32704
 [3] toEPSG(spref::ArchGDAL.ISpatialRef)
   @ ArchGDAL ~/Documents/GitHub/ArchGDAL.jl/src/spatialref.jl:626
 [4] unsafe_convertcrs(#unused#::Type{GeoFormatTypes.EPSG}, crsref::ArchGDAL.ISpatialRef)
   @ ArchGDAL ~/Documents/GitHub/ArchGDAL.jl/src/convert.jl:125
 [5] convert(target::Type{GeoFormatTypes.EPSG}, mode::GeoFormatTypes.CRS, source::GeoFormatTypes.WellKnownText{GeoFormatTypes.CRS})
   @ ArchGDAL ~/Documents/GitHub/ArchGDAL.jl/src/convert.jl:106
 [6] convert(target::Type{GeoFormatTypes.EPSG}, source::GeoFormatTypes.WellKnownText{GeoFormatTypes.CRS}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
   @ GeoFormatTypes ~/.julia/packages/GeoFormatTypes/Lkriv/src/GeoFormatTypes.jl:118
 [7] convert(target::Type{GeoFormatTypes.EPSG}, source::GeoFormatTypes.WellKnownText{GeoFormatTypes.CRS})
   @ GeoFormatTypes ~/.julia/packages/GeoFormatTypes/Lkriv/src/GeoFormatTypes.jl:95
 [8] top-level scope
   @ ~/Documents/GitHub/ArchGDAL.jl/src/Untitled-1.md:10

[alex-s-gardner]

Since the cbinding is over my head, could someone tell be how to properly call 'GDAL.osrgetauthoritycode(spref, epsg)' from within ArchGDAL

[yeesian]

Using the wkt in https://github.com/yeesian/ArchGDAL.jl/issues/390#issue-1805758090 (formatting it for other readers to read the authority codes easily) and following https://gist.github.com/refactor/d41a844c8f5826eca0bceca479a5f0b5, I think it's something like

import GDAL
import ArchGDAL
wkt = """GEOGCS[
    "WGS 84",
    DATUM[
        "WGS_1984",
        SPHEROID[
            "WGS 84",
            6378137,
            298.257223563,
            AUTHORITY["EPSG","7030"]
        ],
        AUTHORITY["EPSG","6326"]
    ],
    PRIMEM[
        "Greenwich",
        0,
        AUTHORITY["EPSG","8901"]
    ],
    UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],
    AXIS["Latitude",NORTH],
    AXIS["Longitude",EAST],
    AUTHORITY["EPSG","4326"]
]
"""
spref = ArchGDAL.importWKT(wkt)
GDAL.osrgetauthoritycode(spref.ptr, "GEOGCS") # returns "4326"
GDAL.osrgetauthoritycode(spref.ptr, "GEOGCS|UNIT") # returns "9122"
GDAL.osrgetauthoritycode(spref.ptr, "GEOGCS|DATUM") # returns "6326"
GDAL.osrgetauthoritycode(spref.ptr, "GEOGCS|DATUM|SPHEROID") # returns "7030"

[visr]

Nice example. One thing to add is that one must not assume that the code an EPSG code, but check this explicitly with OSRGetAuthorityName in GDAL or proj_get_id_auth_name in PROJ.

[alex-s-gardner]

OK, i think I have a working version with this PR

[felixcremer]

Is this done with #391 or is there something missing?

[alex-s-gardner]

Yes, this should be done now.