How to load a raster of a NC file with multiple monthly bands

Specific Domains Geo
,
1

I am using the chelsa_cmip6 script to download climatic data, and the output is a NC file with the data embedded in monthly series, like in this (small) example of precipitations.
How do I load it into a vector of rasters (perhaps better) or a 3 dimensional raster ?
If I load it without parameters using Rasters.Raster I have a strange object without the X dimension:

julia> import NCDatasets;

julia> using Rasters;

julia> ncpath = "chelsa_monthly_precipitations_2014.nc";

julia> ncmap  = Rasters.Raster(ncpath)
╭─────────────────────────────────────────────────╮
│ 2×24 Raster{Union{Missing, Float64},2} lat_bnds │
├─────────────────────────────────────────────────┴─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── dims ┐
  ↓ bnds,
  → Y    Mapped{Float64} [47.004027229800016, 47.01236056310003, …, 47.187360562400016, 47.19569389570003] ForwardOrdered Regular Points
├───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── metadata ┤
  Metadata{Rasters.NCDsource} of Dict{String, Any} with 1 entry:
  "_FillValue" => NaN
├─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── raster ┤
  extent: Extent(bnds = (1, 2), Y = (47.004027229800016, 47.19569389570003))
  missingval: missing
  crs: EPSG:4326
  mappedcrs: EPSG:4326
└───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┘
 47.004   47.0124  47.0207  47.029   47.0374  47.0457  47.054   47.0624  47.0707  47.079   …  47.129   47.1374  47.1457  47.154   47.1624  47.1707  47.179   47.1874  47.1957
 46.0714  46.0797  46.0881  46.0964  46.1047  46.1131  46.1214  46.1297  46.1381  46.1464     46.1964  46.2047  46.2131  46.2214  46.2297  46.2381  46.2464  46.2547  46.2631
 47.9366  47.945   47.9533  47.9616  47.97    47.9783  47.9866  47.995   48.0033  48.0116     48.0616  48.07    48.0783  48.0866  48.095   48.1033  48.1116  48.12    48.1283

julia> ds     = NCDatasets.NCDataset(ncpath)
Dataset: chelsa_monthly_precipitations_2014.nc
Group: /

Dimensions
   lat = 24
   lon = 48
   month = 12
   bnds = 2

Variables
  lat   (24)
    Datatype:    Union{Missing, Float64} (Float64)
    Dimensions:  lat
    Attributes:
     _FillValue           = NaN
     standard_name        = latitude
     long_name            = latitude
     units                = degrees_north

  lon   (48)
    Datatype:    Union{Missing, Float64} (Float64)
    Dimensions:  lon
    Attributes:
     _FillValue           = NaN
     standard_name        = longitude
     long_name            = longitude
     units                = degrees_east

  month   (12)
    Datatype:    Dates.DateTime (Int64)
    Dimensions:  month
    Attributes:
     units                = days since 2014-01-15 00:00:00
     calendar             = proleptic_gregorian

  lat_bnds   (2 × 24)
    Datatype:    Union{Missing, Float64} (Float64)
    Dimensions:  bnds × lat
    Attributes:
     _FillValue           = NaN

  lon_bnds   (2 × 48)
    Datatype:    Union{Missing, Float64} (Float64)
    Dimensions:  bnds × lon
    Attributes:
     _FillValue           = NaN

  pr   (48 × 24 × 12)
    Datatype:    Union{Missing, Float32} (Float32)
    Dimensions:  lon × lat × month
    Attributes:
     _FillValue           = NaN

Global attributes
  coordinates          = lat_bnds lon_bnds


julia> size(ncmap)
(2, 24)

This is how it appears on QGIS (plot based on the first month):

image
image1176×673 161 KB

I can obtain the row data with ds["pr"][:,:,:], perhaps I can use the metadata of the ds object to create a vector of rasters manually using the data of ds["pr"][:,:,month] ??

Python code to get the example data

$ pip install chelsa_cmip6

from chelsa_cmip6.GetClim import chelsa_cmip6
chelsa_cmip6(activity_id='CMIP', 
             table_id='Amon', 
             experiment_id='historical', 
             institution_id='MPI-M', 
             source_id='MPI-ESM1-2-LR', 
             member_id='r1i1p1f1', 
             refps='1981-01-15', 
             refpe='2010-12-15', 
             fefps='2014-01-01', 
             fefpe='2014-12-31', 
             xmin=9.8, 
             xmax=10.2,
             ymin=47.0, 
             ymax=47.2,
             output='~',
             use_esgf=False)
2

Best to make a bug report. But normally you want a RasterStack when there are multiple layers, or specify Raster(path; name=:layerIwant) or you just get the first one. Here the first one found should actually be skipped, it’s a bounds matrix.

1 Like
3

I have tried with “name” but I don’t know what to put on it… (“month” doesn’t seem to work).

I can get a single NCDataSet with, for example, feb = NCDatasets.@select(ds, Dates.month(month) == 2), but then I don’t know how to convert it to a Raster.

Thanks, I’ll open a feature request on the Rasters package.

4

From the structure of your data I think it should be Raster(path; name="pr")which should give you a three dimensional Raster.

5

Thank you… I thought I tried it, instead…