Would be good to add a to_epsg implimentation

[alex-s-gardner]

robustly retrieving the epsg from a crs requires an implementation of pyproj's to_epsg

[visr]

Isn't this already present?

Proj.jl/src/crs.jl

Lines 129 to 133 in d88918c

	function GFT.EPSG(crs::CRS)
	str = proj_get_id_code(crs)
	code = parse(Int, str)
	return GFT.EPSG(code)
	end

Although it looks like this doesn't use proj_get_id_auth_name to check if it is EPSG, so that's a bug.

[alex-s-gardner]

I believe this is a very limited version of to_epsg as it can only find an EPSG code if the Authority and code is contained in the CRS. For example it can't infer the EPSG code from a ProjString or wkt that does not include the authority. This is the same limitation as my ArchGDAL PR. to_epsg has a bunch of logic to try and infer the EPSG and also returns a confidence level that it found the right code.

[alex-s-gardner]

I think we'll need to use Proj.proj_identify

[alex-s-gardner]

Can someone provide an example of how to use Proj.proj_identify, I can't seem to figure out how to assign and retrieve the C objects.

[visr]

This seems to work:

using Proj
import GeoFormatTypes as GFT

# create a CRS from the proj string that represents EPSG:4326
crs = Proj.CRS("+proj=longlat +datum=WGS84 +no_defs +type=crs")
# attempt to identify the corresponding EPSG code using proj_identify
out_confidence = Ref(Ptr{Cint}(C_NULL))
pj_list = Proj.proj_identify(crs, "EPSG", out_confidence)
# TODO use these to check for error and list length
pj_list == C_NULL
Proj.proj_list_get_count(pj_list)
# get the first resulting CRS
pj_1 = Proj.proj_list_get(pj_list, 0)
crs_1 = Proj.CRS(pj_1)
confidence_1 = unsafe_load(out_confidence[])
# use i in unsafe_load(p::Ptr{T}, i::Integer=1) for other indices

# doesn't work: https://github.com/JuliaGeo/Proj.jl/issues/88#issuecomment-1645894753
GFT.EPSG(crs)
# does work:
GFT.EPSG(crs_1) == GFT.EPSG(4326)

• C API docs: https://proj.org/en/9.2/development/reference/functions.html#c.proj_identify
• Proper use of proj_identify in pyproj: https://github.com/pyproj4/pyproj/blob/510c3c2a98cf2f9d1bd14587f865a44245a60330/pyproj/_crs.pyx#L2888

[alex-s-gardner]

Should we introduce this into proj as to_epsg() then use that for GFT.EPSG(crs::CRS)?

[visr]

Yeah perhaps we can create a Proj.identify that is a nicely wrapped version of the C API Proj.proj_identify.

Then for GFT.EPSG(crs::CRS) we can use that function with auth_name set to EPSG and only accepting a confidence of 70 or higher like pyproj.

I just noticed my use of out_confidence was incorrect, that is an array that is filled by proj with the confidence level per entry.

[alex-s-gardner]

In the example above how should out_confidence be initialized then? I don't have a great mental model yet when
calling out to C

[visr]

Hmm yeah I'm also not great with this stuff, but managed to figure it out, see the updated example.

So you need a reference to a pointer of the right type. We can initialize this with a null pointer, PROJ will point it to an array that it allocates itself. After calling the function we can see that out_confidence[] now points somewhere, and use unsafe_load to get the value(s). In this example I got 70. Note the C API docs that explain that we have to free the PROJ allocated array ourselves with Proj.proj_int_list_destroy(out_confidence[]).

[alex-s-gardner]

I'll work on a PR

[alex-s-gardner]

@visr for the higher level identify I think we should just return the top result. I can really see a use case for returning all matches. If that functionality is needed the user can drop a level to proj_identify

[alex-s-gardner]

@visr I think this is good to merge. If you agree I will work on to_EPSG next

[visr]

Left some comments on #92. For to_EPSG you mean GFT.EPSG(crs::CRS) right?