From 6e3bb8b9baf00ecca64bd53d67636f250aa79493 Mon Sep 17 00:00:00 2001 From: Joaquim Date: Mon, 22 Jun 2026 20:17:08 +0100 Subject: [PATCH] Add wkt2epsg and proj2epsg funtions and uncomment some FFT wrapper funs. --- src/GMT.jl | 6 +++-- src/gdal/gdal.jl | 57 ++++++++++++++++++++++++++++++++++++++++++- src/laszip/lazread.jl | 12 ++++----- src/libgmt.jl | 13 +++++----- src/proj_utils.jl | 23 +++++++++++++++++ 5 files changed, 96 insertions(+), 15 deletions(-) diff --git a/src/GMT.jl b/src/GMT.jl index d1c2130f0b..a7ffe9f96f 100644 --- a/src/GMT.jl +++ b/src/GMT.jl @@ -162,7 +162,7 @@ export psimage, psimage!, pslegend, pslegend!, psmask, psmask!, psrose, psrose!, psscale, psscale!, pssolar, pssolar!, psternary, psternary!, pstext, pstext!, pswiggle, pswiggle!, psxy, psxy!, psxyz, psxyz!, regress, resetGMT, rose, rose!, sample1d, scatter, scatter!, scatter3, scatter3!, solar, solar!, analemma, enso, keeling, - sunsetrise, spectrum1d, sphdistance, sphinterpolate, + sunsetrise, spectrum1d, fft1d, sphdistance, sphinterpolate, sphtriangulate, surface, ternary, ternary!, text, text!, text_record, textrepel, annotate, annotate!, trend1d, trend2d, triangulate, gmtsplit, decorated, vector_attrib, wiggle, wiggle!, xyz2grd, gmtbegin, gmtend, gmthelp, subplot, gmtfig, inset, showfig, earthtide, gmt2grd, gravfft, gmtgravmag3d, gravmag3d, grdgravmag3d, gravprisms, grdseamount, parkermag, parkergrav, kovesi, @@ -187,7 +187,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, + geodesic, orthodrome, proj2wkt, setcoords!, setfld!, setcrs!, setsrs!, settimecol!, vecangles, wkt2proj, wkt2epsg, proj2epsg, 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, @@ -336,6 +336,7 @@ include("signalcorr.jl") include("sealand.jl") include("spatial_funs.jl") include("spectrum1d.jl") +include("fft1d.jl") include("sphdistance.jl") include("sphinterpolate.jl") include("sphtriangulate.jl") @@ -420,6 +421,7 @@ using .Laszip gmtwrite(t, [0.0 0; 1 1]) gmtread(t) gmtread(tests("burro_cenora.jpg")) + gmtselect(t, R="0/1/0/1") basemap(R=:d, proj=:guess) coast(R=:g, proj=:guess, W=(level=1, pen=(2, :green)), savefig=tempname()*".ps") rm(t) diff --git a/src/gdal/gdal.jl b/src/gdal/gdal.jl index 5669860375..b831eab1b5 100644 --- a/src/gdal/gdal.jl +++ b/src/gdal/gdal.jl @@ -435,6 +435,52 @@ OSRImportFromEPSG(a1, a2) = acare(ccall((:OSRImportFromEPSG, libgdal), Cint, (pV OSRNewSpatialReference(a1) = acare(ccall((:OSRNewSpatialReference, libgdal), pVoid, (Cstring,), a1)) OSRSetPROJSearchPaths(a1) = acare(ccall((:OSRSetPROJSearchPaths, libgdal), Cvoid, (Ptr{Cstring},), a1)) +""" + OSRAutoIdentifyEPSG(OGRSpatialReferenceH hSRS) -> OGRErr + +Set EPSG authority info if possible. +""" +OSRAutoIdentifyEPSG(a1) = acare(ccall((:OSRAutoIdentifyEPSG, libgdal), Cint, (pVoid,), a1)) + +""" + OSRGetAuthorityName(OGRSpatialReferenceH hSRS, const char * pszTargetKey) -> const char * + +Get the authority name for a node. +""" +function OSRGetAuthorityName(hSRS, pszTargetKey) + acare(ccall((:OSRGetAuthorityName, libgdal), Cstring, (pVoid, Cstring), hSRS, pszTargetKey), false) +end + +""" + OSRGetAuthorityCode(OGRSpatialReferenceH hSRS, const char * pszTargetKey) -> const char * + +Get the authority code for a node. +""" +function OSRGetAuthorityCode(hSRS, pszTargetKey) + acare(ccall((:OSRGetAuthorityCode, libgdal), Cstring, (pVoid, Cstring), hSRS, pszTargetKey), false) +end + + +""" + OSRFindMatches(OGRSpatialReferenceH hSRS, char **papszOptions, int *pnEntries, + int **ppanMatchConfidence) -> OGRSpatialReferenceH * + +Try to identify a match between the passed SRS and a related SRS in a catalog. + +### Parameters +* **hSRS**: SRS to match +* **papszOptions**: NULL terminated list of options or NULL +* **pnEntries**: Output parameter. Number of values in the returned array. +* **ppanMatchConfidence**: Output parameter (or NULL). *ppanMatchConfidence will be allocated to an array of *pnEntries whose values between 0 and 100 indicate the confidence in the match. 100 is the highest confidence level. The array must be freed with CPLFree(). + +### Returns +an array of SRS that match the passed SRS, or NULL. Must be freed with OSRFreeSRSArray() +""" +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 + + function OSRSetAxisMappingStrategy(hSRS, strategy) #(Gdal.GDALVERSION[] < v"3.0.0") && return # This breakes precompile if called from one PrecompileTools call acare(ccall((:OSRSetAxisMappingStrategy, libgdal), Cvoid, (pVoid, UInt32), hSRS, strategy)) @@ -1886,6 +1932,15 @@ end return unsafe_string(projptr[]) 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) + end + + toKML(geom::AbstractGeometry, altitudemode=C_NULL) = OGR_G_ExportToKML(geom.ptr, altitudemode) function importWKT!(spref::AbstractSpatialRef, wktstr::AbstractString) @@ -2653,7 +2708,7 @@ end # --------------------------------- export - bwareaopen, getband, getdriver, getlayer, getproj, getgeom, getgeotransform, toPROJ4, toWKT, importPROJ4, + bwareaopen, getband, getdriver, getlayer, getproj, getgeom, getgeotransform, toPROJ4, toWKT, toEPSG, importPROJ4, importWKT, importEPSG, gdalinfo, gdalwarp, gdaldem, gdaltranslate, gdalgrid, ogr2ogr, gdalvectortranslate, gdalrasterize, gdalbuildvrt, readraster, setgeotransform!, setproj!, destroy, arcellipse, arccircle, delaunay, dither, buffer, centroid, intersection, intersects, polyunion, fromWKT, fillnodata!, fillnodata, diff --git a/src/laszip/lazread.jl b/src/laszip/lazread.jl index 115c226137..dbec7b2ad5 100644 --- a/src/laszip/lazread.jl +++ b/src/laszip/lazread.jl @@ -107,7 +107,7 @@ function lazread(fname::String, out::String, type::DataType, class::Int, startst end xyz = zeros(fType, totalNP, n_col) (occursin('i', argout)) && (intens = zeros(UInt16, totalNP, 1)) - (occursin('c', argout)) && (class = zeros(Int8, totalNP, 1)) + (occursin('c', argout)) && (_class = zeros(UInt8, totalNP, 1)) (occursin('n', argout)) && (n_ret = zeros(Int8, totalNP, 1)) if (n_upper > 0) RGB = zeros(UInt16, totalNP, n_upper) @@ -179,9 +179,9 @@ function lazread(fname::String, out::String, type::DataType, class::Int, startst xyz[k,1] = pt.X; xyz[k,2] = pt.Y; xyz[k,3] = pt.Z intens[k] = pt.intensity if (header.point_data_format > 5) - class[k] = (pt.extended_classification != 0) ? pt.extended_classification : pt.classification + _class[k] = (pt.extended_classification != 0) ? pt.extended_classification : pt.classification else - class[k] = pt.classification + _class[k] = pt.classification end end elseif (argout == "xyzc") @@ -190,9 +190,9 @@ function lazread(fname::String, out::String, type::DataType, class::Int, startst pt = unsafe_load(point[]) xyz[k,1] = pt.X; xyz[k,2] = pt.Y; xyz[k,3] = pt.Z if (header.point_data_format > 5) - class[k] = (pt.extended_classification != 0) ? pt.extended_classification : pt.classification + _class[k] = (pt.extended_classification != 0) ? pt.extended_classification : pt.classification else - class[k] = pt.classification + _class[k] = pt.classification end end elseif ((startswith(argout, "xyz") || startswith(argout, "xy")) && occursin(r"[RGBI]", argout)) @@ -249,7 +249,7 @@ function lazread(fname::String, out::String, type::DataType, class::Int, startst elseif (argout == "xyzi") lasout_types(dsv = [mat2ds(xyz), mat2ds(intens)]) elseif (argout == "xyzc") - lasout_types(dsv = [mat2ds(xyz), mat2ds(class)]) + lasout_types(dsv = [mat2ds(xyz), mat2ds(_class)]) elseif (argout == "xyzti") lasout_types(dsv = [mat2ds(xyz), mat2ds(intens)]) elseif ((startswith(argout, "xyz") || startswith(argout, "xy")) && occursin(r"[RGBI]", argout)) diff --git a/src/libgmt.jl b/src/libgmt.jl index dc9452222f..20d96bfc11 100644 --- a/src/libgmt.jl +++ b/src/libgmt.jl @@ -190,6 +190,11 @@ function GMT_Parse_Common(API::Ptr{Cvoid}, given_options::Ptr{UInt8}, options::P ccall( (:GMT_Parse_Common, libgmt), Cint, (Ptr{Cvoid}, Ptr{UInt8}, Ptr{GMT_OPTION}), API, given_options, options) end +function GMT_Report(API, vlevel::Integer, txt) + ccall((:GMT_Report, libgmt), Cvoid, (Cstring, Cint, Ptr{UInt8}), API, vlevel, txt) +end +=# + function GMT_FFT_Option(API::Ptr{Cvoid}, option::UInt8, dim::UInt32, string::Ptr{UInt8}) ccall( (:GMT_FFT_Option, libgmt), UInt32, (Ptr{Cvoid}, UInt8, UInt32, Ptr{UInt8}), API, option, dim, string) end @@ -208,17 +213,13 @@ end function GMT_FFT_Destroy(API::Ptr{Cvoid}, K::Ptr{Cvoid}) ccall( (:GMT_FFT_Destroy, libgmt), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), API, K) end -function GMT_FFT_1D(API::Ptr{Cvoid}, data::Ptr{Cfloat}, n::UInt16, direction::Cint, mode::UInt32) - ccall( (:GMT_FFT_1D, libgmt), Cint, (Ptr{Cvoid}, Ptr{Cfloat}, UInt16, Cint, UInt32), API, data, n, direction, mode) +function GMT_FFT_1D(API::Ptr{Cvoid}, data::Ptr{Cfloat}, n::Integer, direction::Cint, mode::UInt32) + ccall( (:GMT_FFT_1D, libgmt), Cint, (Ptr{Cvoid}, Ptr{Cfloat}, UInt64, Cint, UInt32), API, data, UInt64(n), direction, mode) end function GMT_FFT_2D(API::Ptr{Cvoid}, data::Ptr{Cfloat}, nx::UInt32, ny::UInt32, direction::Cint, mode::UInt32) ccall( (:GMT_FFT_2D, libgmt), Cint, (Ptr{Cvoid}, Ptr{Cfloat}, UInt32, UInt32, Cint, UInt32), API, data, nx, ny, direction, mode) end -function GMT_Report(API, vlevel::Integer, txt) - ccall((:GMT_Report, libgmt), Cvoid, (Cstring, Cint, Ptr{UInt8}), API, vlevel, txt) -end =# - function GMT_Encode_Options(V_API::Ptr{Cvoid}, _module, n_argin::Int, head::Ref{Ptr{GMT_OPTION}}, n) ccall((:GMT_Encode_Options, libgmt), Ptr{GMT_RESOURCE}, (Cstring, Ptr{UInt8}, Int32, Ref{Ptr{GMT_OPTION}}, Ptr{UInt32}), V_API, _module, n_argin, head, n) diff --git a/src/proj_utils.jl b/src/proj_utils.jl index 18224a8f6a..c5483d3d48 100644 --- a/src/proj_utils.jl +++ b/src/proj_utils.jl @@ -792,6 +792,7 @@ function proj_ellipsoid_get_parameters(ellipsoid::Ptr{Cvoid}=NULL, proj_ptr::Ptr semi_major[], inv_flattening[], semi_minor[] end +# ----------------------------------------------------------------------------------------- """ proj2wkt(proj4_str::String, pretty::Bool=false) @@ -802,6 +803,26 @@ function proj2wkt(proj4_str::String; pretty::Bool=false) toWKT(importPROJ4(proj4_str), pretty) end +# ----------------------------------------------------------------------------------------- +""" + proj2epsg(proj4_str::String) + +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))) +end + +# ----------------------------------------------------------------------------------------- +""" + wkt2epsg(wkt_str::String) + +Convert a WKT SRS string into the corresponding EPSG code (if it exists). +""" +wkt2epsg(wkt_str::String) = toEPSG(importWKT(wkt_str)) + +# ----------------------------------------------------------------------------------------- """ wkt2proj(wkt_str::String) @@ -809,6 +830,7 @@ Convert a WKT SRS string into the PROJ4 form. """ wkt2proj(wkt_str::String) = toPROJ4(importWKT(wkt_str)) +# ----------------------------------------------------------------------------------------- """ epsg2proj(code::Integer) @@ -816,6 +838,7 @@ Convert a EPSG code into the PROJ4 form. """ epsg2proj(code::Integer) = toPROJ4(importEPSG(code)) +# ----------------------------------------------------------------------------------------- """ epsg2wkt(code::Integer, pretty::Bool=false)