From b6862c86bc7441f8080c8356db1f985728d6098a Mon Sep 17 00:00:00 2001 From: Joaquim Date: Mon, 19 Oct 2026 16:48:43 +0100 Subject: [PATCH] Add a wrapper to gdalsrsinfo and use it to convert proj2epsg --- src/GMT.jl | 2 +- src/gdal/gdal.jl | 40 ++++++++++++++++++++++++++++++++++++---- src/proj_utils.jl | 25 ++++++++++++++++++++++++- test/test_proj4.jl | 16 ++++++++++++++++ 4 files changed, 77 insertions(+), 6 deletions(-) diff --git a/src/GMT.jl b/src/GMT.jl index 684cbe230d..d17713ab64 100644 --- a/src/GMT.jl +++ b/src/GMT.jl @@ -194,7 +194,7 @@ export wkbMultiLineString, wkbMultiPolygon, wkbGeometryCollection, wkbPoint25D, wkbLineString25D, wkbPolygon25D, wkbMultiPoint25D, wkbMultiLineString25D, wkbMultiPolygon25D, wkbGeometryCollection25D, bezier, buffergeo, circgeo, epsg2proj, epsg2wkt, geod, invgeod, loxodrome, loxodrome_direct, loxodrome_inverse, montage, - geodesic, orthodrome, proj2wkt, setcoords!, setfld!, setcrs!, setsrs!, settimecol!, vecangles, wkt2proj, wkt2epsg, proj2epsg, + geodesic, orthodrome, proj2wkt, setcoords!, setfld!, setcrs!, setsrs!, settimecol!, vecangles, wkt2proj, wkt2epsg, proj2epsg, srsmatches, inbbox, randgeo, colorzones!, rasterzones!, rasterzones, lelandshade, texture_img, crop, doy2date, date2doy, choropleth, fourcolors, getdcw, ISOtime2unix, median, mean, quantile, std, nanmean, nanstd, skipnan, zonal_statistics, zonal_stats, autocor, autocor!, autocov, autocov!, conv, yeardecimal, xcorr, xcov, add2PSfile, append2fig, isoutlier, linearfitxy, diff --git a/src/gdal/gdal.jl b/src/gdal/gdal.jl index dc12b351e0..e3a6037bca 100644 --- a/src/gdal/gdal.jl +++ b/src/gdal/gdal.jl @@ -486,6 +486,7 @@ an array of SRS that match the passed SRS, or NULL. Must be freed with OSRFreeSR function OSRFindMatches(hSRS, papszOptions, pnEntries, ppanMatchConfidence) acare(ccall((:OSRFindMatches, libgdal), Ptr{pVoid}, (pVoid, Ptr{Cstring}, Ptr{Cint}, Ptr{Ptr{Cint}}), hSRS, papszOptions, pnEntries, ppanMatchConfidence)) end +OSRFreeSRSArray(a1) = acare(ccall((:OSRFreeSRSArray, libgdal), Cvoid, (Ptr{pVoid},), a1)) function OSRSetAxisMappingStrategy(hSRS, strategy) @@ -1945,10 +1946,41 @@ end function toEPSG(spref::AbstractSpatialRef)::Int (spref.ptr == C_NULL) && return 0 - (OSRAutoIdentifyEPSG(spref.ptr) != OGRERR_NONE) && error("Failed to identify EPSG code for this SRS") - result = OSRGetAuthorityCode(spref.ptr, C_NULL) - (result === "") && error("Failed to get EPSG code") - return parse(Int, result) + if (OSRAutoIdentifyEPSG(spref.ptr) == OGRERR_NONE) + result = OSRGetAuthorityCode(spref.ptr, C_NULL) + (result !== "") && return parse(Int, result) + end + # AutoIdentifyEPSG only reads authority nodes already in the SRS. A bare PROJ4 (or a WKT made + # from one) has none, so search the PROJ database instead, as gdalsrsinfo -o epsg does. + for (code, _) in findmatches(spref) + startswith(code, "EPSG:") && return parse(Int, code[6:end]) + end + error("Failed to identify EPSG code for this SRS") + end + + """ + findmatches(spref::AbstractSpatialRef) -> Vector{Tuple{String,Int}} + + Search the PROJ database for the registered systems matching `spref` (OSRFindMatches). Returns + `("AUTHORITY:CODE", confidence)` pairs, best first, confidence in 0-100. Empty if nothing matches. + """ + function findmatches(spref::AbstractSpatialRef)::Vector{Tuple{String,Int}} + out = Tuple{String,Int}[] + (spref.ptr == C_NULL) && return out + n, conf = Ref{Cint}(0), Ref{Ptr{Cint}}(C_NULL) + arr = OSRFindMatches(spref.ptr, C_NULL, n, conf) + try + for k = 1:n[] + h = unsafe_load(arr, k) + auth, code = OSRGetAuthorityName(h, C_NULL), OSRGetAuthorityCode(h, C_NULL) + (auth === "" || code === "") && continue + push!(out, (auth * ":" * code, Int(unsafe_load(conf[], k)))) + end + finally + (arr != C_NULL) && OSRFreeSRSArray(arr) + (conf[] != C_NULL) && VSIFree(conf[]) + end + return out end diff --git a/src/proj_utils.jl b/src/proj_utils.jl index c5483d3d48..bfd97e8ac8 100644 --- a/src/proj_utils.jl +++ b/src/proj_utils.jl @@ -811,7 +811,7 @@ Convert a PROJ4 string into the corresponding EPSG code (if it exists). """ function proj2epsg(proj4_str::String) !startswith(proj4_str, "+proj=") && error("$(proj4_str) is not a valid proj4 string") - toEPSG(importWKT(proj2wkt(proj4_str))) + toEPSG(importPROJ4(proj4_str)) # Not via a WKT round-trip: that loses the EPSG:4326 match of a plain longlat end # ----------------------------------------------------------------------------------------- @@ -822,6 +822,29 @@ Convert a WKT SRS string into the corresponding EPSG code (if it exists). """ wkt2epsg(wkt_str::String) = toEPSG(importWKT(wkt_str)) +# ----------------------------------------------------------------------------------------- +""" + srsmatches(srs::String) -> Vector{Tuple{String,Int}} + +Find the registered reference systems that match `srs` (a PROJ4 string, a WKT string or "EPSG:n"), +as `gdalsrsinfo -e` does. Returns `("AUTHORITY:CODE", confidence)` pairs, best first, with the +confidence in 0-100. Matches can come from other authorities than EPSG (OGC, IGNF, ESRI, ...). +Returns an empty vector when nothing matches (e.g. a bare `+proj=moll`). + +### Example +```julia +srsmatches("+proj=utm +zone=29 +datum=WGS84") +1-element Vector{Tuple{String, Int64}}: + ("EPSG:32629", 70) +``` +""" +function srsmatches(srs::String) + s = strip(srs) + sr = startswith(s, "+proj=") ? importPROJ4(s) : + startswith(uppercase(s), "EPSG:") ? importEPSG(parse(Int, s[6:end])) : importWKT(s) + Gdal.findmatches(sr) +end + # ----------------------------------------------------------------------------------------- """ wkt2proj(wkt_str::String) diff --git a/test/test_proj4.jl b/test/test_proj4.jl index c118424465..57f5bc0051 100644 --- a/test/test_proj4.jl +++ b/test/test_proj4.jl @@ -27,6 +27,22 @@ prj = epsg2proj(4326) proj2wkt(prj) wkt2proj(wkt) + + # EPSG from a bare PROJ4/WKT (no authority nodes): the OSRFindMatches fallback in toEPSG + @test proj2epsg("+proj=utm +zone=29 +datum=WGS84") == 32629 + @test proj2epsg("+proj=utm +zone=29 +south +datum=WGS84") == 32729 + @test proj2epsg("+proj=merc") == 3395 + @test proj2epsg("+proj=longlat +datum=WGS84 +no_defs") == 4326 + @test wkt2epsg(proj2wkt("+proj=utm +zone=29 +datum=WGS84")) == 32629 + @test wkt2epsg(epsg2wkt(32629)) == 32629 + @test_throws ErrorException proj2epsg("+proj=moll +lon_0=0") + @test_throws ErrorException proj2epsg("lixo") + @test srsmatches("+proj=utm +zone=29 +datum=WGS84") == [("EPSG:32629", 70)] + @test ("EPSG:4326", 70) in srsmatches("+proj=longlat +datum=WGS84 +no_defs") + @test srsmatches("EPSG:32629") == [("EPSG:32629", 100)] + @test srsmatches("epsg:3763")[1] == ("EPSG:3763", 100) + @test srsmatches(epsg2wkt(32629))[1] == ("EPSG:32629", 100) + @test isempty(srsmatches("+proj=moll +lon_0=0")) loxo = loxodrome_direct(0,0,45, 10000) loxo = loxodrome([0 0; 30 50], step=500, unit=:k); loxo = loxodrome(mat2ds([0 0; 30 50]), step=500, unit=:k);