Add a wrapper to gdalsrsinfo and use it to convert proj2epsg

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,

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
 
 

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)

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);