Add multi file support for Grib and NetCDF files

[aasdelat]

I am aware that Rasters.jl uses GRIBDatasets.jl and NCDatasets.jl in order to deal with Gribs and NetCDF files, respectively.

Recently, the ability to make a single dataset from a list of files, has been added to these two libraries (see this thread).

Is it possible that Rasters.jl takes advantage of this new feature in order to make a Raster or RasterStack (but not a RasterSeries) from a bunch of files?

[lazarusA]

In yax we have: https://juliadatacubes.github.io/YAXArrays.jl/stable/UserGuide/read.html#open-mfdataset , but I haven't tested it with Grib files.

[rafaqz]

What does it do exactly?

We have always been able to do merge(RasterStack(filea), RasterStack(fileb)) to get a single RasterStack. How would this be different?

Ok maybe I get it a bit more: its specifically concatenation? In Rasters you can just do cat(stack1, stack2; dims=Ti) and it will lazily concatenate if the stacks are lazy. Is that what you want?

[rafaqz]

Also how is it different from combine(rasterseries) ? I clearly don't get it...

[aasdelat]

I will use an example:
I have two grib files with the same structure. They have two dimensions: "values", that are integer indices of spatial points, and Ti, that are timestamps. The timestamps of the second file are continuation of the timestamps of the first file. And both files have three variables: longitude, latitude and msl (mean sea level pressure).
I can read both files as a RasterSeries:
dataset = RasterSeries(input_files, :valid_time, child=RasterStack, lazy=true)
Here are structures of both files:

julia> dataset[1]
╭────────────────────────╮
│ 338880×744 RasterStack │
├────────────────────────┴─────────────────────────────────────────────────────────────────────────── dims ┐
  ↓ values,
  → Ti     Sampled{DateTime} [1980-12-01T00:00:00, …, 1980-12-31T23:00:00] ForwardOrdered Irregular Points
├────────────────────────────────────────────────────────────────────────────────────────────────── layers ┤
  :longitude eltype: Union{Missing, Float64} dims: values size: 338880
  :latitude  eltype: Union{Missing, Float64} dims: values size: 338880
  :msl       eltype: Union{Missing, Float64} dims: values, Ti size: 338880×744
├──────────────────────────────────────────────────────────────────────────────────────────────────────────┴─────────────────────────── metadata ┐
  Metadata{Rasters.GRIBsource} of Dict{String, Any} with 6 entries:
  "edition"           => "1"
  "source"            => "/home/antonio/Documentos/datos/ERA5/ERA5_N320_1h_sl_sfc_NH15S_1978_2023_msl/ERA5_N320_1h_sl_sfc_NH15S_1980-12-01_1980-1…
  "centreDescription" => "European Centre for Medium-Range Weather Forecasts"
  "centre"            => "ecmf"
  "subCentre"         => "0"
  "Conventions"       => "CF-1.7"
├──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── raster ┤
  extent: Extent(values = (1, 338880), Ti = (DateTime("1980-12-01T00:00:00"), DateTime("1980-12-31T23:00:00")))
  missingval: missing
  filename: 
└────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘

julia> dataset[2]
╭────────────────────────╮
│ 338880×744 RasterStack │
├────────────────────────┴─────────────────────────────────────────────────────────────────────────── dims ┐
  ↓ values,
  → Ti     Sampled{DateTime} [1981-01-01T00:00:00, …, 1981-01-31T23:00:00] ForwardOrdered Irregular Points
├────────────────────────────────────────────────────────────────────────────────────────────────── layers ┤
  :longitude eltype: Union{Missing, Float64} dims: values size: 338880
  :latitude  eltype: Union{Missing, Float64} dims: values size: 338880
  :msl       eltype: Union{Missing, Float64} dims: values, Ti size: 338880×744
├──────────────────────────────────────────────────────────────────────────────────────────────────────────┴─────────────────────────── metadata ┐
  Metadata{Rasters.GRIBsource} of Dict{String, Any} with 6 entries:
  "edition"           => "1"
  "source"            => "/home/antonio/Documentos/datos/ERA5/ERA5_N320_1h_sl_sfc_NH15S_1978_2023_msl/ERA5_N320_1h_sl_sfc_NH15S_1981-01-01_1981-0…
  "centreDescription" => "European Centre for Medium-Range Weather Forecasts"
  "centre"            => "ecmf"
  "subCentre"         => "0"
  "Conventions"       => "CF-1.7"
├──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── raster ┤
  extent: Extent(values = (1, 338880), Ti = (DateTime("1981-01-01T00:00:00"), DateTime("1981-01-31T23:00:00")))
  missingval: missing
  filename: 
└────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘

Note that the time dimension of the second file is the continuation of the one of the first file.
So, I want to have only one RasterStack object that results from the concatenation of both files.
If I use combine, I get:

julia> combine(dataset[1:2], dims=:Ti)
╭────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────╮
│ 2-element DimArray{RasterStack{(:longitude, :latitude, :msl), @NamedTuple{longitude::Union{Missing, Float64}, latitude::Union{Missing, Float64}, msl::Union{Missing, Float64}}, 2, @NamedTuple{longitude::DiskArrays.ConcatDiskArray{Union{Missing, Float64}, 1, Vector{Raster{Union{Missing, Float64}, 1, Tuple{Dim{:values, DimensionalData.Dimensions.Lookups.NoLookup{Base.OneTo{Int64}}}}, Tuple{}, Rasters.FileArray{Rasters.GRIBsource, Union{Missing, Float64}, 1, Symbol, Nothing, DiskArrays.GridChunks{1, Tuple{DiskArrays.RegularChunks}}, DiskArrays.Unchunked}, Symbol, DimensionalData.Dimensions.Lookups.Metadata{Rasters.GRIBsource, Dict{String, Any}}, Missing}}}, latitude::DiskArrays.ConcatDiskArray{Union{Missing, Float64}, 1, Vector{Raster{Union{Missing, Float64}, 1, Tuple{Dim{:values, DimensionalData.Dimensions.Lookups.NoLookup{Base.OneTo{Int64}}}}, Tuple{}, Rasters.FileArray{Rasters.GRIBsource, Union{Missing, Float64}, 1, Symbol, Nothing, DiskArrays.GridChunks{1, Tuple{DiskArrays.RegularChunks}}, DiskArrays.Unchunked}, Symbol, DimensionalData.Dimensions.Lookups.Metadata{Rasters.GRIBsource, Dict{String, Any}}, Missing}}}, msl::DiskArrays.ConcatDiskArray{Union{Missing, Float64}, 2, Matrix{Raster{Union{Missing, Float64}, 2, Tuple{Dim{:values, DimensionalData.Dimensions.Lookups.NoLookup{Base.OneTo{Int64}}}, Ti{DimensionalData.Dimensions.Lookups.Sampled{DateTime, Vector{DateTime}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Irregular{Tuple{Nothing, Nothing}}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.Metadata{Rasters.GRIBsource, Dict{String, Any}}}}}, Tuple{}, Rasters.FileArray{Rasters.GRIBsource, Union{Missing, Float64}, 2, Symbol, Nothing, DiskArrays.GridChunks{2, Tuple{DiskArrays.RegularChunks, DiskArrays.RegularChunks}}, DiskArrays.Unchunked}, Symbol, DimensionalData.Dimensions.Lookups.Metadata{Rasters.GRIBsource, Dict{String, Any}}, Missing}}}}, Tuple{Dim{:values, DimensionalData.Dimensions.Lookups.NoLookup{Base.OneTo{Int64}}}, Ti{DimensionalData.Dimensions.Lookups.Sampled{DateTime, Vector{DateTime}, DimensionalData.Dimensions.Lookups.ForwardOrdered, DimensionalData.Dimensions.Lookups.Irregular{Tuple{Nothing, Nothing}}, DimensionalData.Dimensions.Lookups.Points, DimensionalData.Dimensions.Lookups.Metadata{Rasters.GRIBsource, Dict{String, Any}}}}}, Tuple{}, @NamedTuple{longitude::Tuple{Dim{:values, Colon}}, latitude::Tuple{Dim{:values, Colon}}, msl::Tuple{Dim{:values, Colon}, Ti{Colon}}}, DimensionalData.Dimensions.Lookups.Metadata{Rasters.GRIBsource, Dict{String, Any}}, @NamedTuple{longitude::DimensionalData.Dimensions.Lookups.Metadata{Rasters.GRIBsource, Dict{String, Any}}, latitude::DimensionalData.Dimensions.Lookups.Metadata{Rasters.GRIBsource, Dict{String, Any}}, msl::DimensionalData.Dimensions.Lookups.Metadata{Rasters.GRIBsource, Dict{String, Any}}}, Missing},1} │
├────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── dims ┤
  ↓ valid_time
└────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
 RasterStack{(:longitude, :latitude, :msl), @NamedTuple{longitude::Union{Missing, Float64}, latitude::Union{Missing, Float64}, msl::Union{Missing, Float64}}, 2, @NamedTuple{longitude::ConcatDiskArray{Union{Missing, Float64}, 1, Vector{Raster{Union{Missing, Float64}, 1, Tuple{Dim{:values, NoLookup{OneTo{Int64}}}}, Tuple{}, FileArray{GRIBsource, Union{Missing, Float64}, 1, Symbol, Nothing, GridChunks{1, Tuple{RegularChunks}}, Unchunked}, Symbol, Metadata{GRIBsource, Dict{String, Any}}, Missing}}}, latitude::ConcatDiskArray{Union{Missing, Float64}, 1, Vector{Raster{Union{Missing, Float64}, 1, Tuple{Dim{:values, NoLookup{OneTo{Int64}}}}, Tuple{}, FileArray{GRIBsource, Union{Missing, Float64}, 1, Symbol, Nothing, GridChunks{1, Tuple{RegularChunks}}, Unchunked}, Symbol, Metadata{GRIBsource, Dict{String, Any}}, Missing}}}, msl::ConcatDiskArray{Union{Missing, Float64}, 2, Matrix{Raster{Union{Missing, Float64}, 2, Tuple{Dim{:values, NoLookup{OneTo{Int64}}}, Ti{Sampled{DateTime, Vector{DateTime}, ForwardOrdered, Irregular{Tuple{Nothing, Nothing}}, Points, Metadata{GRIBsource, Dict{String, Any}}}}}, Tuple{}, FileArray{GRIBsource, Union{Missing, Float64}, 2, Symbol, Nothing, GridChunks{2, Tuple{RegularChunks, RegularChunks}}, Unchunked}, Symbol, Metadata{GRIBsource, Dict{String, Any}}, Missing}}}}, Tuple{Dim{:values, NoLookup{OneTo{Int64}}}, Ti{Sampled{DateTime, Vector{DateTime}, ForwardOrdered, Irregular{Tuple{Nothing, Nothing}}, Points, Metadata{GRIBsource, Dict{String, Any}}}}}, Tuple{}, @NamedTuple{longitude::Tuple{Dim{:values, Colon}}, latitude::Tuple{Dim{:values, Colon}}, msl::Tuple{Dim{:values, Colon}, Ti{Colon}}}, Metadata{GRIBsource, Dict{String, Any}}, @NamedTuple{longitude::Metadata{GRIBsource, Dict{String, Any}}, latitude::Metadata{GRIBsource, Dict{String, Any}}, msl::Metadata{GRIBsource, Dict{String, Any}}}, Missing}((longitude = 338880-element ConcatDiskArray{Union{Missing, Float64}, 1, Vector{Raster{Union{Missing, Float64}, 1, Tuple{Dim{:values, NoLookup{OneTo{Int64}}}}, Tuple{}, FileArray{GRIBsource, Union{Missing, Float64}, 1, Symbol, Nothing, GridChunks{1, Tuple{RegularChunks}}, Unchunked}, Symbol, Metadata{GRIBsource, Dict{String, Any}}, Missing}}}, latitude = 338880-element ConcatDiskArray{Union{Missing, Float64}, 1, Vector{Raster{Union{Missing, Float64}, 1, Tuple{Dim{:values, NoLookup{OneTo{Int64}}}}, Tuple{}, FileArray{GRIBsource, Union{Missing, Float64}, 1, Symbol, Nothing, GridChunks{1, Tuple{RegularChunks}}, Unchunked}, Symbol, Metadata{GRIBsource, Dict{String, Any}}, Missing}}}, msl = 338880×744 ConcatDiskArray{Union{Missing, Float64}, 2, Matrix{Raster{Union{Missing, Float64}, 2, Tuple{Dim{:values, NoLookup{OneTo{Int64}}}, Ti{Sampled{DateTime, Vector{DateTime}, ForwardOrdered, Irregular{Tuple{Nothing, Nothing}}, Points, Metadata{GRIBsource, Dict{String, Any}}}}}, Tuple{}, FileArray{GRIBsource, Union{Missing, Float64}, 2, Symbol, Nothing, GridChunks{2, Tuple{RegularChunks, RegularChunks}}, Unchunked}, Symbol, Metadata{GRIBsource, Dict{String, Any}}, Missing}}}), ↓ values,
→ Ti     Sampled{DateTime} [1980-12-01T00:00:00, …, 1980-12-31T23:00:00] ForwardOrdered Irregular Points, (), (longitude = ↓ values, latitude = ↓ values, msl = ↓ values, → Ti), Metadata{GRIBsource, Dict{String, Any}}(Dict{String, Any}("edition"=>"1", "source"=>"/home/antonio/Documentos/datos/ERA5/ERA5_N320_1h_sl_sfc_NH15S_1978_2023_msl/ERA5_N320_1h_sl_sfc_NH15S_1980-12-01_1980-12-31_msl.grib", "centreDescription"=>"European Centre for Medium-Range Weather Forecasts", "centre"=>"ecmf", "subCentre"=>"0", "Conventions"=>"CF-1.7")), (longitude = Metadata{GRIBsource, Dict{String, Any}}(Dict{String, Any}("units"=>"degrees_east", "long_name"=>"longitude", "standard_name"=>"longitude")), latitude = Metadata{GRIBsource, Dict{String, Any}}(Dict{String, Any}("units"=>"degrees_north", "long_name"=>"latitude", "standard_name"=>"latitude")), msl = Metadata{GRIBsource, Dict{String, Any}}(Dict{String, Any}("units"=>"Pa", "long_name"=>"Mean sea level pressure", "standard_name"=>"air_pressure_at_mean_sea_level"))), missing)
 RasterStack{(:longitude, :latitude, :msl), @NamedTuple{longitude::Union{Missing, Float64}, latitude::Union{Missing, Float64}, msl::Union{Missing, Float64}}, 2, @NamedTuple{longitude::ConcatDiskArray{Union{Missing, Float64}, 1, Vector{Raster{Union{Missing, Float64}, 1, Tuple{Dim{:values, NoLookup{OneTo{Int64}}}}, Tuple{}, FileArray{GRIBsource, Union{Missing, Float64}, 1, Symbol, Nothing, GridChunks{1, Tuple{RegularChunks}}, Unchunked}, Symbol, Metadata{GRIBsource, Dict{String, Any}}, Missing}}}, latitude::ConcatDiskArray{Union{Missing, Float64}, 1, Vector{Raster{Union{Missing, Float64}, 1, Tuple{Dim{:values, NoLookup{OneTo{Int64}}}}, Tuple{}, FileArray{GRIBsource, Union{Missing, Float64}, 1, Symbol, Nothing, GridChunks{1, Tuple{RegularChunks}}, Unchunked}, Symbol, Metadata{GRIBsource, Dict{String, Any}}, Missing}}}, msl::ConcatDiskArray{Union{Missing, Float64}, 2, Matrix{Raster{Union{Missing, Float64}, 2, Tuple{Dim{:values, NoLookup{OneTo{Int64}}}, Ti{Sampled{DateTime, Vector{DateTime}, ForwardOrdered, Irregular{Tuple{Nothing, Nothing}}, Points, Metadata{GRIBsource, Dict{String, Any}}}}}, Tuple{}, FileArray{GRIBsource, Union{Missing, Float64}, 2, Symbol, Nothing, GridChunks{2, Tuple{RegularChunks, RegularChunks}}, Unchunked}, Symbol, Metadata{GRIBsource, Dict{String, Any}}, Missing}}}}, Tuple{Dim{:values, NoLookup{OneTo{Int64}}}, Ti{Sampled{DateTime, Vector{DateTime}, ForwardOrdered, Irregular{Tuple{Nothing, Nothing}}, Points, Metadata{GRIBsource, Dict{String, Any}}}}}, Tuple{}, @NamedTuple{longitude::Tuple{Dim{:values, Colon}}, latitude::Tuple{Dim{:values, Colon}}, msl::Tuple{Dim{:values, Colon}, Ti{Colon}}}, Metadata{GRIBsource, Dict{String, Any}}, @NamedTuple{longitude::Metadata{GRIBsource, Dict{String, Any}}, latitude::Metadata{GRIBsource, Dict{String, Any}}, msl::Metadata{GRIBsource, Dict{String, Any}}}, Missing}((longitude = 338880-element ConcatDiskArray{Union{Missing, Float64}, 1, Vector{Raster{Union{Missing, Float64}, 1, Tuple{Dim{:values, NoLookup{OneTo{Int64}}}}, Tuple{}, FileArray{GRIBsource, Union{Missing, Float64}, 1, Symbol, Nothing, GridChunks{1, Tuple{RegularChunks}}, Unchunked}, Symbol, Metadata{GRIBsource, Dict{String, Any}}, Missing}}}, latitude = 338880-element ConcatDiskArray{Union{Missing, Float64}, 1, Vector{Raster{Union{Missing, Float64}, 1, Tuple{Dim{:values, NoLookup{OneTo{Int64}}}}, Tuple{}, FileArray{GRIBsource, Union{Missing, Float64}, 1, Symbol, Nothing, GridChunks{1, Tuple{RegularChunks}}, Unchunked}, Symbol, Metadata{GRIBsource, Dict{String, Any}}, Missing}}}, msl = 338880×744 ConcatDiskArray{Union{Missing, Float64}, 2, Matrix{Raster{Union{Missing, Float64}, 2, Tuple{Dim{:values, NoLookup{OneTo{Int64}}}, Ti{Sampled{DateTime, Vector{DateTime}, ForwardOrdered, Irregular{Tuple{Nothing, Nothing}}, Points, Metadata{GRIBsource, Dict{String, Any}}}}}, Tuple{}, FileArray{GRIBsource, Union{Missing, Float64}, 2, Symbol, Nothing, GridChunks{2, Tuple{RegularChunks, RegularChunks}}, Unchunked}, Symbol, Metadata{GRIBsource, Dict{String, Any}}, Missing}}}), ↓ values,
→ Ti     Sampled{DateTime} [1981-01-01T00:00:00, …, 1981-01-31T23:00:00] ForwardOrdered Irregular Points, (), (longitude = ↓ values, latitude = ↓ values, msl = ↓ values, → Ti), Metadata{GRIBsource, Dict{String, Any}}(Dict{String, Any}("edition"=>"1", "source"=>"/home/antonio/Documentos/datos/ERA5/ERA5_N320_1h_sl_sfc_NH15S_1978_2023_msl/ERA5_N320_1h_sl_sfc_NH15S_1981-01-01_1981-01-31_msl.grib", "centreDescription"=>"European Centre for Medium-Range Weather Forecasts", "centre"=>"ecmf", "subCentre"=>"0", "Conventions"=>"CF-1.7")), (longitude = Metadata{GRIBsource, Dict{String, Any}}(Dict{String, Any}("units"=>"degrees_east", "long_name"=>"longitude", "standard_name"=>"longitude")), latitude = Metadata{GRIBsource, Dict{String, Any}}(Dict{String, Any}("units"=>"degrees_north", "long_name"=>"latitude", "standard_name"=>"latitude")), msl = Metadata{GRIBsource, Dict{String, Any}}(Dict{String, Any}("units"=>"Pa", "long_name"=>"Mean sea level pressure", "standard_name"=>"air_pressure_at_mean_sea_level"))), missing)

And I do not understand what this is.

On the other hand, I can use cat. But, logically, cat will concatenate only the variable with a time dimension: msl, but neither longitude nor latitude. So, first, I have to select the variable that will be concatenated along the time asis (msl). Second, I have to pass to cat, the new (time)axis values, id est, I have to concatenate, first, the time axes of the two files:

time_name = :Ti
var_name = :msl
lon_name = :longitude
lat_name = :latitude

new_time_dim = cat(
    [ lookup(dataset[i], time_name) for i in 1:length(dataset)]...,
    dims=1
    )#cat
new_ds_msl = cat(
    [ dataset[i][var_name] for i in 1:length(dataset)]..., # We have selected the :msl variable
    dims=Ti(new_time_dim)
    )#cat
# We have now the msl variable (or Raster), in new_ds_msl (a Raster object), but I have to
# reincorporate longitude and latitude, resulting a RasterStack object:
dataset = RasterStack(
    ( dataset[1][lon_name], dataset[1][lat_name], new_ds_var ),
    name = (lon_name, lat_name, var_name)
    )#RasterStack

And I have what I need:

julia> dataset
╭─────────────────────────╮
│ 338880×1488 RasterStack │
├─────────────────────────┴────────────────────────────────────────────────────────────────────────── dims ┐
  ↓ values,
  → Ti     Sampled{DateTime} [1980-12-01T00:00:00, …, 1981-01-31T23:00:00] ForwardOrdered Irregular Points
├────────────────────────────────────────────────────────────────────────────────────────────────── layers ┤
  :longitude eltype: Union{Missing, Float64} dims: values size: 338880
  :latitude  eltype: Union{Missing, Float64} dims: values size: 338880
  :msl       eltype: Union{Missing, Float64} dims: values, Ti size: 338880×1488
├────────────────────────────────────────────────────────────────────────────────────────────────── raster ┤
  extent: Extent(values = (1, 338880), Ti = (DateTime("1980-12-01T00:00:00"), DateTime("1981-01-31T23:00:00")))
  missingval: missing
  filename: 
└──────────────────────────────────────────────────────────────────────────────────────────────────────────┘

And I think that you will tell me: "You have what you want, so What's the problem then?"
The problem is that this all staff is a bit cumbersome and should be simplified in a simple function like xarray.open_mfdataset. This makes the code more readable and easier to be maintained.

So, although I can do what I want with the current implementation, I think it would be very advisable to have a xarray.open_mfdataset like function to make things simpler and more intuitive.

[rafaqz]

Huh.

It seems the reason this doesn't work easily is because your data (lat, lon) is really a dimension, not data. This needs to load as a Vector data cube with a point lookup. (We can do this but it's not yet being detected on loading the data)

Then the points would be the same for both stacks, and cat would just work on the time dimension as I suggested.

You could also use spatial selectors, which you currently can't do.

(@asinghvi17 here's a simple vector cube example in the wild)

Do you have a link to the data?

(And there is no way I would say that much code is ok! This should all be one or two lines. The question is which lines... I don't want to encourage dimensions as data like you have, youre losing plotting and half the tools in the package with that setup)

[rafaqz]

@aasdelat a link to that file would really help identify how we can automatically load it as a dimension.

I'm not sure if its following CF point standards, or if it is custom data that is just organised like that.

[rafaqz]

Also of note is :

#917

We can just accept already loaded datasets from CDM.jl and turn them into RasterStacks

[rafaqz]

FYI you can probably just do Raster(multifile) and RasterStack(multifile) now that those constructors on AbstractDataset are merged.

[aasdelat]

Sorry form my (so long) delay.
I will answer in order:

The files are 0,5 Gb each and are retrieved from the Copernicus Data Store. You have to be registered in order to download data from here (and also you have to know what data you want and how to retrieve them), but I share you the two files by my Mega account:
https://mega.nz/file/4gZE3A4K#xWFXWlHil2lAM5NpNAC-WTSEIa6PhYK8BTokOoW7iwE
https://mega.nz/file/k9Zl0bCA#RFMs-F7rgjc1EG8MWnEL2yjmJoqR6o7ITNzk_hL_T7M

Sorry, but the link you provide in

(@asinghvi17 here's a simple vector cube example in the wild)

is to the profile of Anshul Singhvi , not to an example of a vector cube.

The organization of the spatial structure information in these files is very common in the atmospheric sciences, although some people prefer to interpolate them into regular lat/lon grids, but I prefer to work with the "pristine" data. This spatial organization is not very friendly for the plotting tools out there. Panply manages them very well, but lacks of configuring options. Metview seems to do calculations with them, but cannot plot them properly. It would be great to plot them with GeoMakie (I have not tested this, due to my lack of knowledge about it).

RasterStack(mfds) or Raster(mfds) does not work:

files = ["ERA5_N320_1h_sl_sfc_NH15S_1980-12-01_1980-12-31_msl.grib", "ERA5_N320_1h_sl_sfc_NH15S_1981-01-01_1981-01-31_msl.grib"]
ds = GRIBDataset.(files)
mfds = MFDataset(ds,aggdim = "valid_time")
rmfds = Raster(mfds)
ERROR: MethodError: no method matching sourcetrait(::MFDataset{GRIBDataset{Float64, 2, Missing}, 1, String})

Closest candidates are:
  sourcetrait(::Symbol)
   @ Rasters ~/.julia/packages/Rasters/mJsGi/src/sources/sources.jl:100
  sourcetrait(::Type{<:Rasters.Source})
   @ Rasters ~/.julia/packages/Rasters/mJsGi/src/sources/sources.jl:99
  sourcetrait(::CommonDataModel.CFVariable)
   @ Rasters ~/.julia/packages/Rasters/mJsGi/src/sources/commondatamodel.jl:39
  ...

Stacktrace:
 [1] _raster(ds::MFDataset{GRIBDataset{Float64, 2, Missing}, 1, String})
   @ Rasters ~/.julia/packages/Rasters/mJsGi/src/array.jl:334
 [2] Raster(ds::MFDataset{GRIBDataset{Float64, 2, Missing}, 1, String}; kw::@Kwargs{})
   @ Rasters ~/.julia/packages/Rasters/mJsGi/src/array.jl:333
 [3] top-level scope
   @ REPL[39]:1

OR

rmfds = RasterStack(mfds)
ERROR: MethodError: no method matching sourcetrait(::MFDataset{GRIBDataset{Float64, 2, Missing}, 1, String})

Closest candidates are:
  sourcetrait(::Symbol)
   @ Rasters ~/.julia/packages/Rasters/mJsGi/src/sources/sources.jl:100
  sourcetrait(::Type{<:Rasters.Source})
   @ Rasters ~/.julia/packages/Rasters/mJsGi/src/sources/sources.jl:99
  sourcetrait(::CommonDataModel.CFVariable)
   @ Rasters ~/.julia/packages/Rasters/mJsGi/src/sources/commondatamodel.jl:39
  ...

Stacktrace:
 [1] _cdmlookup(ds::MFDataset{…}, dimname::String, D::Type, crs::Rasters.NoKW, mappedcrs::Rasters.NoKW)
   @ Rasters ~/.julia/packages/Rasters/mJsGi/src/sources/commondatamodel.jl:277
 [2] _cdmdim(ds::MFDataset{GRIBDataset{…}, 1, String}, dimname::String, crs::Rasters.NoKW, mappedcrs::Rasters.NoKW)
   @ Rasters ~/.julia/packages/Rasters/mJsGi/src/sources/commondatamodel.jl:226
 [3] _dimdict(ds::MFDataset{GRIBDataset{Float64, 2, Missing}, 1, String}, crs::Rasters.NoKW, mappedcrs::Rasters.NoKW)
   @ Rasters ~/.julia/packages/Rasters/mJsGi/src/sources/commondatamodel.jl:173
 [4] RasterStack(ds::MFDataset{…}; filename::String, source::Rasters.NoKW, dims::Rasters.NoKW, refdims::Tuple{}, name::Rasters.NoKW, group::Rasters.NoKW, metadata::Rasters.NoKW, layermetadata::Rasters.NoKW, layerdims::Rasters.NoKW, missingval::Rasters.NoKW, crs::Rasters.NoKW, mappedcrs::Rasters.NoKW, coerce::typeof(convert), checkmem::Bool, scaled::Rasters.NoKW, lazy::Bool, verbose::Bool, raw::Bool, kw::@Kwargs{})
   @ Rasters ~/.julia/packages/Rasters/mJsGi/src/stack.jl:476
 [5] RasterStack(ds::MFDataset{GRIBDataset{Float64, 2, Missing}, 1, String})
   @ Rasters ~/.julia/packages/Rasters/mJsGi/src/stack.jl:451
 [6] top-level scope
   @ REPL[38]:1
Some type information was truncated. Use `show(err)` to see complete types.

[rafaqz]

Thats actually not the worst error to see, we simply need to get the trait for GRIB from the inner dataset. I'm guessing it cant hold mixed dataset types?

[asinghvi17]

Looking at the struct definition it seems like there's nothing explicitly against it but we can always error, if they are not the same sourcetrait... can get e.g. only(unique(sourcetrait, inner_datasets)).

@aasdelat do you have an image of what you want to be plotted from this dataset? Plotting as points is simple but there are many ways to plot these kinds of datasets on a map, maybe you want to interpolate?

Also, what exactly do you mean by performing math on these cubes? If it's the exact same cube, with the exact same lat and long values, then it makes sense - but otherwise, what does it mean to perform math on long/lat vars that don't match up?

[asinghvi17]

also, @rafaqz this is what the attributes of the dataset look like:

Dataset: /Users/anshul/Downloads/ERA5_N320_1h_sl_sfc_NH15S_1980-12-01_1980-12-31_msl.grib
Group: /

Dimensions
   values = 338880
   valid_time = 744

Variables
  valid_time   (744)
    Datatype:    Dates.DateTime (Int64)
    Dimensions:  valid_time
    Attributes:
     units                = seconds since 1970-01-01T00:00:00
     calendar             = proleptic_gregorian
     long_name            = time
     standard_name        = time

  longitude   (338880)
    Datatype:    Float64 (Float64)
    Dimensions:  values
    Attributes:
     units                = degrees_east
     long_name            = longitude
     standard_name        = longitude

  latitude   (338880)
    Datatype:    Float64 (Float64)
    Dimensions:  values
    Attributes:
     units                = degrees_north
     long_name            = latitude
     standard_name        = latitude

  msl   (338880 × 744)
    Datatype:    Union{Missing, Float64} (Float64)
    Dimensions:  values × valid_time
    Attributes:
     units                = Pa
     long_name            = Mean sea level pressure
     standard_name        = air_pressure_at_mean_sea_level

Global attributes
  edition              = 1
  source               = /Users/anshul/Downloads/ERA5_N320_1h_sl_sfc_NH15S_1980-12-01_1980-12-31_msl.grib
  centreDescription    = European Centre for Medium-Range Weather Forecasts
  centre               = ecmf
  subCentre            = 0
  Conventions          = CF-1.7

so there's no way for us to know what the dimensions are there, this doesn't seem to conform to CF standards either as there is no "coordinates" attribute to msl.

[rafaqz]

Yeah it can still work without coordinates, we know Lat and Lon are connected to values and are standard Lattitude and longitude, while there is no values variable

[rafaqz]

So it's a Point Geometry lookup

(And I guess a 3d xy scatter over time)