From e50580673fd03bab752102f9f68399ad426dc562 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Thu, 10 Apr 2025 16:17:57 -0400 Subject: [PATCH 1/2] load CDM CRS'es if they are available in string form --- src/sources/commondatamodel.jl | 46 ++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) diff --git a/src/sources/commondatamodel.jl b/src/sources/commondatamodel.jl index f0fab3101..59b06c5d5 100644 --- a/src/sources/commondatamodel.jl +++ b/src/sources/commondatamodel.jl @@ -168,6 +168,10 @@ function _layers(ds::AbstractDataset, names, group) end function _dims(var::AbstractVariable{<:Any,N}, crs=nokw, mappedcrs=nokw) where N + if isnokw(crs) + # Get CRS by hook or by crook - start at the variable and work our way up + crs = _get_crs_by_hook_or_crook(CDM.dataset(var), var, CDM.attribs(var)) + end dimnames = CDM.dimnames(var) ntuple(Val(N)) do i _cdmdim(CDM.dataset(var), dimnames[i], crs, mappedcrs) @@ -339,6 +343,48 @@ function _cdmlookup( end return _cdmlookup(D, index, order, span, sampling, metadata, crs, mappedcrs) end + +function _get_crs_by_hook_or_crook(ds::AbstractDataset, var, attr) + # Get CRS by hook or by crook - start at the variable and work our way up + crs = nokw + if haskey(attr, "grid_mapping") + if haskey(ds, attr["grid_mapping"]) + crs = _crs_from_cf_attr(CDM.attribs(ds[attr["grid_mapping"]])) + !isnokw(crs) && return crs + else + crs = _crs_from_cf_attr(attr) + !isnokw(crs) && return crs + end + elseif haskey(CDM.attribs(var), "grid_mapping") + var_attr = CDM.attribs(var) + if haskey(ds, var_attr["grid_mapping"]) + crs = _crs_from_cf_attr(CDM.attribs(ds[var_attr["grid_mapping"]])) + !isnokw(crs) && return crs + else + crs = _crs_from_cf_attr(var_attr) + !isnokw(crs) && return crs + end + else# nothing locally...check global attributes as a last resort + return _crs_from_cf_attr(attr) + end + +end + +function _crs_from_cf_attr(attr) + # we can't parse CF crs fully, BUT + # many CF implementations will include either `crs_wkt`, `spatial_epsg`, or `proj4string` + # which we do know the types of + if haskey(attr, "crs_wkt") + return ESRIWellKnownText(GeoFormatTypes.CRS(), attr["crs_wkt"]) + elseif haskey(attr, "spatial_epsg") + return EPSG(parse(Int, attr["spatial_epsg"])) + elseif haskey(attr, "proj4string") + return ProjString(attr["proj4string"]) + end + return nokw # in the worst case return nothing +end + + # For X and Y use a Mapped <: AbstractSampled lookup function _cdmlookup( D::Type{<:Union{<:XDim,<:YDim}}, index, order::Order, span, sampling, metadata, crs, mappedcrs From 02bec44ddfe19c86a045edec4f619e214ff78113 Mon Sep 17 00:00:00 2001 From: Anshul Singhvi Date: Mon, 14 Apr 2025 13:33:48 -0400 Subject: [PATCH 2/2] return WKT2, not ESRIWKT because WKT2 is guaranteed to be backwards compatible with WKT1 --- src/sources/commondatamodel.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/sources/commondatamodel.jl b/src/sources/commondatamodel.jl index 59b06c5d5..68ab7a705 100644 --- a/src/sources/commondatamodel.jl +++ b/src/sources/commondatamodel.jl @@ -375,7 +375,7 @@ function _crs_from_cf_attr(attr) # many CF implementations will include either `crs_wkt`, `spatial_epsg`, or `proj4string` # which we do know the types of if haskey(attr, "crs_wkt") - return ESRIWellKnownText(GeoFormatTypes.CRS(), attr["crs_wkt"]) + return WellKnownText2(GeoFormatTypes.CRS(), attr["crs_wkt"]) elseif haskey(attr, "spatial_epsg") return EPSG(parse(Int, attr["spatial_epsg"])) elseif haskey(attr, "proj4string")