Proj.identify (#92)

* add high level interface to proj_identify * add tests for Proj.identify * move identify from libproj to crs * address comments on Proj.identify * update tests * free pj_list * final improvements - standardize docstring - avoid intermediate Vectors - free out_confidence * run JuliaFormatter --------- Co-authored-by: Martijn Visser <[email protected]>
JuliaGeo · Aug 30, 2023 · 5beab68 · 5beab68
1 parent f67081e
commit 5beab68
Show file tree

Hide file tree

Showing 3 changed files with 62 additions and 8 deletions.
diff --git a/src/coord.jl b/src/coord.jl
@@ -88,12 +88,19 @@ end
 function Transformation(
     source_crs::GFT.CoordinateReferenceSystemFormat,
     target_crs::GFT.CoordinateReferenceSystemFormat;
-    always_xy::Bool=false,
-    direction::PJ_DIRECTION=PJ_FWD,
-    area::Ptr{PJ_AREA}=C_NULL,
-    ctx::Ptr{PJ_CONTEXT}=C_NULL
+    always_xy::Bool = false,
+    direction::PJ_DIRECTION = PJ_FWD,
+    area::Ptr{PJ_AREA} = C_NULL,
+    ctx::Ptr{PJ_CONTEXT} = C_NULL,
 )
-    return Transformation(CRS(source_crs).pj, CRS(target_crs).pj; always_xy, direction, area, ctx)
+    return Transformation(
+        CRS(source_crs).pj,
+        CRS(target_crs).pj;
+        always_xy,
+        direction,
+        area,
+        ctx,
+    )
 end
 
 function Transformation(

diff --git a/src/crs.jl b/src/crs.jl
@@ -1,7 +1,7 @@
 """
     CRS(crs)
 
-Create a CRS. `crs` can be:
+Create a coordinate reference system. `crs` can be:
 - a proj-string,
 - a WKT string,
 - an object code (like "EPSG:4326", "urn:ogc:def:crs:EPSG::4326", "urn:ogc:def:coordinateOperation:EPSG::1671"),
@@ -137,3 +137,39 @@ Base.convert(::Type{CRS}, crs::GFT.CoordinateReferenceSystemFormat) = CRS(crs)
 
 # Maybe enable later, based on https://github.com/JuliaGeo/GeoFormatTypes.jl/issues/21
 # Base.convert(T::Type{<:GFT.CoordinateReferenceSystemFormat}, crs::GFT.CoordinateReferenceSystemFormat) = T(CRS(crs))
+
+"""
+    identify(crs::CRS; auth_name = nothing)
+
+Returns a list of matching reference CRS and confidence values (0-100).
+
+# Arguments
+- `crs::CRS`: Coordinate reference system
+- `auth_name=nothing`: Authority name, or nothing for all authorities (e.g. "EPSG")
+"""
+function identify(
+    crs::CRS;
+    auth_name = nothing,
+)::Vector{@NamedTuple{crs::CRS, confidence::Int32}}
+
+    out_confidence = Ref(Ptr{Cint}(C_NULL))
+    if isnothing(auth_name)
+        # set authority to C_NULL
+        auth_name = C_NULL
+    end
+
+    pj_list = proj_identify(crs, auth_name, out_confidence)
+    list = NamedTuple{(:crs, :confidence),Tuple{CRS,Int32}}[]
+
+    # was a match found?
+    if pj_list != C_NULL
+        n = proj_list_get_count(pj_list)
+        for i = 1:n
+            crs = CRS(proj_list_get(pj_list, i - 1))
+            confidence = unsafe_load(out_confidence[], i)
+            push!(list, (; crs, confidence))
+        end
+        proj_int_list_destroy(out_confidence[])
+    end
+    return list
+end
diff --git a/test/libproj.jl b/test/libproj.jl
@@ -202,15 +202,15 @@ end
     source_crs = GFT.EPSG("EPSG:4326")
     target_crs = GFT.EPSG("EPSG:32628")
 
-    trans = Proj.Transformation(source_crs, target_crs, always_xy=false)
+    trans = Proj.Transformation(source_crs, target_crs, always_xy = false)
     info = Proj.proj_pj_info(trans.pj)
     description1 = unsafe_string(info.definition)
 
     #crs as txt
     source_crs = "EPSG:4326"
     target_crs = "EPSG:32628"
 
-    trans = Proj.Transformation(source_crs, target_crs, always_xy=false)
+    trans = Proj.Transformation(source_crs, target_crs, always_xy = false)
     info = Proj.proj_pj_info(trans.pj)
     description2 = unsafe_string(info.definition)
 
@@ -477,3 +477,14 @@ end
     # Maybe enable later, based on https://github.com/JuliaGeo/GeoFormatTypes.jl/issues/21
     # @test convert(GFT.ProjString, gftcrs) == GFT.ProjString("+proj=longlat +datum=WGS84 +no_defs +type=crs")
 end
+
+@testset "Proj.identify" begin
+    crs = Proj.CRS("+proj=longlat +datum=WGS84 +no_defs +type=crs")
+    identity = Proj.identify(crs)
+    @test identity[1].confidence == 70
+    @test length(identity) == 4
+
+    identity = Proj.identify(crs; auth_name = "EPSG")
+    @test GFT.EPSG(identity[1].crs) == GFT.EPSG(4326)
+    @test identity[1].confidence == 70
+end