JuliaDataCubes / YAXArrays.jl

Yet Another XArray-like Julia package
https://juliadatacubes.github.io/YAXArrays.jl/
Other
101 stars 18 forks source link

Opening rotated gridded dataset fails to load grid mapping definition #329

Closed Balinus closed 1 year ago

Balinus commented 1 year ago

Hello!

I am trying to play with a rotated grid dataset. YAXArrays is able to load the dataset, but there is missing the grid mapping definition and also misses to load latitude and longitude 2D matrices. This makes it impossible to use the mapping definition for further reprojecting of the data or creating maps with the lat-long matrices.

Here's a MWE.

Cheers!

# using NCDatasets
# using URIs
using YAXArrays
using NetCDF
# using Zarr

# URL RDRS on Pavics
url = "https://pavics.ouranos.ca/twitcher/ows/proxy/thredds/dodsC/datasets/reanalyses/day_RDRSv2.1_NAM.ncml"

# Lecture de la structure
ds = Cube(url)

# Returns 
 706×800×14245×13 YAXArray{Float32,4} with dimensions: 
 Dim{:rlon} Sampled{Float32} 324.60278f0:0.09000002f0:388.0528f0 ForwardOrdered Regular Points,
 Dim{:rlat} Sampled{Float32} Float32[-46.170002, -46.08, …, 25.65, 25.74] ForwardOrdered Irregular Points,
 Ti Sampled{DateTime} DateTime[1980-01-01T00:00:00, …, 2018-12-31T00:00:00] ForwardOrdered Irregular Points,
 Dim{:Variable} Categorical{String} String[huss, hurs, …, tasmin, tdps] Unordered
 units: K
Total size: 389.64 GB

ds.Variable

# Returns
Dim{:Variable} Categorical{String} Unordered
wrapping: 13-element Vector{String}:
 "huss"
"hurs"
"tasmax"
"pr"
"psl"
"tdpsmin"
"hursmax"
"tdpsmax"
"ps"
"hursmin"
"tas"
"tasmin"
"tdps"

ds.axes
# Returns
Dim{:rlon} Sampled{Float32} 324.60278f0:0.09000002f0:388.0528f0 ForwardOrdered Regular Points,
Dim{:rlat} Sampled{Float32} Float32[-46.170002, -46.08, …, 25.65, 25.74] ForwardOrdered Irregular Points,
Ti Sampled{DateTime} DateTime[1980-01-01T00:00:00, …, 2018-12-31T00:00:00] ForwardOrdered Irregular Points,
Dim{:Variable} Categorical{String} String[huss, hurs, …, tasmin, tdps] Unordered

Compared to xarray's output (lat, lon and rotated_pole are available in Coordinates)

import threddsclient
import xarray as xr
import numpy as np

url = 'https://pavics.ouranos.ca/twitcher/ows/proxy/thredds/catalog/datasets/reanalyses/catalog.html'

datasets = [ds for ds in threddsclient.crawl(url) if 'day_RDRSv2.1_NAM.ncml' in ds.name]

assert datasets[0].name
ds = xr.open_dataset(datasets[0].opendap_url())

# Returns 
xarray.Dataset
Dimensions: rlat: 800 rlon: 706 time: 14245
Coordinates:
rlat (rlat) float32 -46.17 -46.08 ... 25.65 25.74
rlon (rlon) float32 324.6 324.7 324.8 ... 388.0 388.1
rotated_pole () float32 ...
time (time) datetime64[ns] 1980-01-01 ... 2018-12-31
lat (rlat, rlon) float32 ...
lon (rlat, rlon) float32 ...

Data variables:
tas (time, rlat, rlon) float32 ...
tasmin (time, rlat, rlon) float32 ...
tasmax (time, rlat, rlon) float32 ...
pr (time, rlat, rlon) float32 ...
tdpsmin (time, rlat, rlon) float32 ...
tdpsmax (time, rlat, rlon) float32 ...
tdps (time, rlat, rlon) float32 ...
psl (time, rlat, rlon) float32 ...
ps (time, rlat, rlon) float32 ...
huss (time, rlat, rlon) float32 ...
hursmin (time, rlat, rlon) float32 ...
hursmax (time, rlat, rlon) float32 ...
hurs (time, rlat, rlon) float32 ...
Balinus commented 1 year ago

Humm, a simple error on my side. I should have used open_dataset(url) and not Cube(url).

Thanks!