Add support for loading CRS in CDM formats

src/extra/cdm.jl

@@ -19,10 +19,9 @@
   cnames = filter(!isnothing, [xname, yname, zname])
   isempty(cnames) && error("coordinates not found")
   coords = map(nm -> _var2array(ds[nm], lazy), cnames)
+  N = length(cnames)
 
-  grid = if all(ndims(a) == 1 for a in coords)
-    RectilinearGrid(coords...)
-  else
+  if !all(ndims(a) == 1 for a in coords)
     error("invalid grid arrays")
   end
 
@@ -35,6 +34,38 @@
     issetequal(vdims, cnames) || issetequal(vdims, dnames)
   end
 
+  # get grid mapping (CRS)
+  gridmappings = map(vnames) do name
+    var = ds[name]
+    attribs = CDM.attribnames(var)
+    if "grid_mapping" ∈ attribs
+      CDM.attrib(var, "grid_mapping")
+    else
+      nothing
+    end
+  end
+
+  # convert grid mapping to CRS
+  crs = if all(isnothing, gridmappings)
+    nothing
+  else
+    if !allequal(gridmappings)
+      error("all variables must have the same CRS")
+    end
+
+    gridmapping = first(gridmappings)
+    _gm2crs(ds[gridmapping])
+  end
+
+  # construct grid with CRS and Manifold
+  C = isnothing(crs) ? Cartesian{NoDatum,N,Met{Float64}} : crs
+  grid = if C <: LatLon
+    lons, lats = coords
+    RectilinearGrid{🌐,C}(lats, lons)
+  else
+    RectilinearGrid{𝔼{N},C}(coords...)
+  end
+
   vtable = if isempty(vnames)
     nothing
   else
@@ -107,6 +138,18 @@
   end
 end
 
+# reference of ellipsoid names: https://raw.githubusercontent.com/wiki/cf-convention/cf-conventions/csv/ellipsoid.csv
+const ELLIP2DATUM = Dict(
+  "WGS 84" => WGS84Latest,
+  "GRS 1980" => ITRFLatest,
+  "Airy 1830" => OSGB36,
+  "Airy Modified 1849" => Ire65,
+  "Bessel 1841" => Hermannskogel,
+  "International 1924" => NZGD1949,
+  "Clarke 1880 (IGN)" => Carthage,
+  "GRS 1967 Modified" => SAD69
+)
+
 function _gribdataset(fname)
   if Sys.iswindows()
     error("loading GRIB files is currently not supported on Windows")
@@ -169,3 +212,49 @@
     name
   end
 end
+
+function _gm2crs(gridmapping)
+  attribs = CDM.attribnames(gridmapping)
+
+  # get datum from reference ellipsoid
+  D = if "reference_ellipsoid_name" ∈ attribs
+    ellip = CDM.attrib(gridmapping, "reference_ellipsoid_name")
+    ELLIP2DATUM[ellip]
+  else
+    WGS84Latest
+  end
+
+  # shift parameters
+  function shift()
+    lonₒ = "longitude_of_central_meridian" ∈ attribs ? CDM.attrib(gridmapping, "longitude_of_central_meridian") : 0.0
+    xₒ = "false_easting" ∈ attribs ? CDM.attrib(gridmapping, "false_easting") : 0.0
+    yₒ = "false_northing" ∈ attribs ? CDM.attrib(gridmapping, "false_northing") : 0.0
+    CoordRefSystems.Shift(lonₒ=lonₒ * u"°", xₒ=xₒ * u"m", yₒ=yₒ * u"m")
+  end
+
+  # parse CRS type and parameters
+  # reference: https://cfconventions.org/cf-conventions/cf-conventions.html#appendix-grid-mappings
+  gmname = CDM.attrib(gridmapping, "grid_mapping_name")
+  if gmname == "latitude_longitude"
+    LatLon{D,Deg{Float64}}
+  elseif gmname == "lambert_cylindrical_equal_area"
+    latₜₛ = if "standard_parallel" ∈ attribs
+      CDM.attrib(gridmapping, "standard_parallel")
+    elseif "scale_factor_at_projection_origin" ∈ attribs
+      CDM.attrib(gridmapping, "scale_factor_at_projection_origin")
+    end
+    CoordRefSystems.EqualAreaCylindrical{latₜₛ * u"°",D,shift(),Met{Float64}}
+  elseif gmname == "mercator"
+    Mercator{D,shift(),Met{Float64}}
+  elseif gmname == "orthographic"
+    latₒ = CDM.attrib(gridmapping, "latitude_of_projection_origin")
+    Mode = CoordRefSystems.EllipticalMode
+    CoordRefSystems.Orthographic{Mode,latₒ * u"°",D,shift(),Met{Float64}}
+  elseif gmname == "transverse_mercator"
+    k₀ = CDM.attrib(gridmapping, "scale_factor_at_central_meridian")
+    latₒ = CDM.attrib(gridmapping, "latitude_of_projection_origin")
+    TransverseMercator{k₀,latₒ * u"°",D,shift(),Met{Float64}}
+  else
+    nothing
+  end
+end

src/utils.jl

@@ -2,6 +2,10 @@
 # Licensed under the MIT License. See LICENSE in the project root.
 # ------------------------------------------------------------------
 
+# helper type alias
+const Met{T} = Quantity{T,u"𝐋",typeof(u"m")}
+const Deg{T} = Quantity{T,NoDims,typeof(u"°")}
+
 function asgeotable(table)
   crs = GI.crs(table)
   cols = Tables.columns(table)

test/io/netcdf.jl

@@ -30,6 +30,25 @@
     @test length(first(vtable.tempanomaly)) == 100
     @test length(vtable.data) == nvertices(gtb.geometry)
     @test eltype(vtable.data) <: Float64
+
+    # CRS
+    file = joinpath(datadir, "test_latlon.nc")
+    gtb = GeoIO.load(file)
+    @test gtb.geometry isa RectilinearGrid
+    @test crs(gtb.geometry) <: LatLon
+    @test datum(crs(gtb.geometry)) === WGS84Latest
+
+    file = joinpath(datadir, "test_latlon_itrf.nc")
+    gtb = GeoIO.load(file)
+    @test gtb.geometry isa RectilinearGrid
+    @test crs(gtb.geometry) <: LatLon
+    @test datum(crs(gtb.geometry)) === ITRFLatest
+
+    file = joinpath(datadir, "test_utm_north_32.nc")
+    gtb = GeoIO.load(file)
+    @test gtb.geometry isa RectilinearGrid
+    @test crs(gtb.geometry) <: TransverseMercator
+    @test datum(crs(gtb.geometry)) === WGS84Latest
   end
 
   @testset "save" begin