Converting a temporal DatasetCollection to NetCDF¶

A DatasetCollection is a time-series of aligned rasters — conceptually a (time, band, y, x) cube backed by a folder of GeoTIFFs. Writing it to a single NetCDF collapses that folder into one self-describing CF file with a real time dimension, which is far easier to move around and open in xarray, QGIS, or Panoply than a pile of dated .tif files.

Input	Output
folder of dated GeoTIFFs	one NetCDF with a time dimension of length T

The whole trip is DatasetCollection.read_multiple_files(folder, with_order=True).to_netcdf(path) — with_order parses the date in each file name so the cube gets a real calendar axis (details below). This notebook covers the multi-timestep case. For single-file raster⇄NetCDF conversions (one GeoTIFF's bands to variables and back), see Converting between GeoTIFF and NetCDF.

Setup¶

In [1]:

%matplotlib inline

import tempfile
from pathlib import Path

import pandas as pd

DATA = Path('../../../examples/data')
WORK = Path(
    tempfile.mkdtemp(prefix='pyramids-collection-nc-')
)  # scratch dir for outputs
DATA.is_dir(), WORK.is_dir()

%matplotlib inline import tempfile from pathlib import Path import pandas as pd DATA = Path('../../../examples/data') WORK = Path( tempfile.mkdtemp(prefix='pyramids-collection-nc-') ) # scratch dir for outputs DATA.is_dir(), WORK.is_dir()

Out[1]:

(True, True)

In [2]:

from pyramids.dataset import DatasetCollection
from pyramids.netcdf import NetCDF
from pyramids.netcdf.plot_options import Selectors

from pyramids.dataset import DatasetCollection from pyramids.netcdf import NetCDF from pyramids.netcdf.plot_options import Selectors

Build a temporal collection¶

read_multiple_files reads a folder of rasters into one stack. With with_order=True it parses the date in each file name (here MSWEP_1979.01.01.tif → file_name_data_fmt='%Y.%m.%d') and orders the timesteps by it, so the collection gets a real calendar time axis rather than an arbitrary file order. The sample folder holds six daily MSWEP rainfall rasters.

In [3]:

folder = DATA / 'geotiff' / 'raster-folder'
sorted(p.name for p in folder.glob('*.tif'))

folder = DATA / 'geotiff' / 'raster-folder' sorted(p.name for p in folder.glob('*.tif'))

Out[3]:

['MSWEP_1979.01.01.tif',
 'MSWEP_1979.01.02.tif',
 'MSWEP_1979.01.03.tif',
 'MSWEP_1979.01.04.tif',
 'MSWEP_1979.01.05.tif',
 'MSWEP_1979.01.06.tif']

In [4]:

col = DatasetCollection.read_multiple_files(
    str(folder), with_order=True, file_name_data_fmt='%Y.%m.%d'
)
col.time_length, col.meta.band_names, col.meta.shape

col = DatasetCollection.read_multiple_files( str(folder), with_order=True, file_name_data_fmt='%Y.%m.%d' ) col.time_length, col.meta.band_names, col.meta.shape

/tmp/ipykernel_2690/1555281888.py:1: DeprecationWarning: DatasetCollection.read_multiple_files is deprecated; use from_files(path, glob=..., date_format=...) instead.
  col = DatasetCollection.read_multiple_files(

Out[4]:

(6, ('Band_1',), (1, 125, 93))

The parsed dates become the collection's time axis — this is what to_netcdf writes out as the NetCDF time coordinate.

In [5]:

times = list(col.time)
times

times = list(col.time) times

Out[5]:

[datetime.datetime(1979, 1, 1, 0, 0),
 datetime.datetime(1979, 1, 2, 0, 0),
 datetime.datetime(1979, 1, 3, 0, 0),
 datetime.datetime(1979, 1, 4, 0, 0),
 datetime.datetime(1979, 1, 5, 0, 0),
 datetime.datetime(1979, 1, 6, 0, 0)]

Look at the stack¶

iloc(i) pulls timestep i out as a plain Dataset, so we can plot it. Compare the first and last days to see the rainfall field change over the week.

In [6]:

col.iloc(0).plot(band=0, title=f'Timestep 0 — {times[0]:%Y-%m-%d}')

col.iloc(0).plot(band=0, title=f'Timestep 0 — {times[0]:%Y-%m-%d}')

Out[6]:

<cleopatra.glyphs.gridded.array_glyph.ArrayGlyph at 0x7f655f258ad0>

In [7]:

col.iloc(col.time_length - 1).plot(
    band=0, title=f'Timestep {col.time_length - 1} — {times[-1]:%Y-%m-%d}'
)

col.iloc(col.time_length - 1).plot( band=0, title=f'Timestep {col.time_length - 1} — {times[-1]:%Y-%m-%d}' )

Out[7]:

<cleopatra.glyphs.gridded.array_glyph.ArrayGlyph at 0x7f65570dc190>

Convert to a single NetCDF¶

One call. The calendar time axis is written automatically — you do not pass it. With the default var_per_band=True, each band becomes its own CF variable and the band names carry through (this single-band stack yields one variable, Band_1).

In [8]:

out = WORK / 'mswep_cube.nc'
col.to_netcdf(out)
out.stat().st_size  # bytes on disk

out = WORK / 'mswep_cube.nc' col.to_netcdf(out) out.stat().st_size # bytes on disk

Out[8]:

299629

Read it back¶

NetCDF.read_file reopens the file through GDAL (no netcdf4 / h5netcdf engine needed). The time dimension and its length survive the round-trip, and each variable is shaped (time, y, x).

In [9]:

nc = NetCDF.read_file(out)
nc.variable_names, nc.epsg

nc = NetCDF.read_file(out) nc.variable_names, nc.epsg

Out[9]:

(['Band_1'], 4647)

In [10]:

nc.dimension_names, nc.dimension_sizes

nc.dimension_names, nc.dimension_sizes

Out[10]:

(['time', 'y', 'x'], {'time': 6, 'y': 125, 'x': 93})

In [11]:

nc.get_variable('Band_1').shape  # (time, y, x)

nc.get_variable('Band_1').shape # (time, y, x)

Out[11]:

(6, 125, 93)

The calendar axis round-trips too: time_stamp decodes the stored CF time coordinate back to the original dates.

In [12]:

nc.time_stamp

nc.time_stamp

/home/runner/work/pyramids/pyramids/.pixi/envs/docs/lib/python3.14/site-packages/osgeo/gdal.py:13234: RuntimeWarning: dimension #0 (time) is not a Time or Vertical dimension.
  return _gdal.Open(*args)

Out[12]:

['1979-01-01',
 '1979-01-02',
 '1979-01-03',
 '1979-01-04',
 '1979-01-05',
 '1979-01-06']

Plot a timestep straight from the reopened NetCDF to confirm the pixels match the source. NetCDF plotting selects along named dimensions, so pick the time slice with Selectors(isel={'time': i}). The source's no-data value round-trips as well, so the fill is masked automatically — no exclude_value needed.

In [13]:

nc.plot(
    variable='Band_1',
    selectors=Selectors(isel={'time': 0}),
    title='Reopened NetCDF — timestep 0',
)

nc.plot( variable='Band_1', selectors=Selectors(isel={'time': 0}), title='Reopened NetCDF — timestep 0', )

Out[13]:

<cleopatra.glyphs.gridded.array_glyph.ArrayGlyph at 0x7f6556fc4e10>

Give the cube your own time axis¶

When the file names carry no parseable date — or you simply want a different calendar — pass time_coords explicitly. Its length must equal time_length.

In [14]:

dated = WORK / 'mswep_2020.nc'
col.to_netcdf(
    dated, time_coords=pd.date_range('2020-01-01', periods=col.time_length, freq='D')
)
NetCDF.read_file(dated).time_stamp

dated = WORK / 'mswep_2020.nc' col.to_netcdf( dated, time_coords=pd.date_range('2020-01-01', periods=col.time_length, freq='D') ) NetCDF.read_file(dated).time_stamp

/home/runner/work/pyramids/pyramids/.pixi/envs/docs/lib/python3.14/site-packages/osgeo/gdal.py:13234: RuntimeWarning: dimension #0 (time) is not a Time or Vertical dimension.
  return _gdal.Open(*args)

Out[14]:

['2020-01-01',
 '2020-01-02',
 '2020-01-03',
 '2020-01-04',
 '2020-01-05',
 '2020-01-06']

One 4-D variable instead of one-per-band¶

For cubes with many bands, var_per_band=False writes a single data variable carrying an extra band dimension, instead of one variable per band.

In [15]:

cube4d = WORK / 'mswep_4d.nc'
col.to_netcdf(cube4d, var_per_band=False)
nc4 = NetCDF.read_file(cube4d)
nc4.variable_names, nc4.dimension_sizes

cube4d = WORK / 'mswep_4d.nc' col.to_netcdf(cube4d, var_per_band=False) nc4 = NetCDF.read_file(cube4d) nc4.variable_names, nc4.dimension_sizes

Out[15]:

(['data'], {'time': 6, 'band': 1, 'y': 125, 'x': 93})

Notes¶

• to_netcdf needs xarray installed — it assembles an xarray.Dataset in memory, then writes it through GDAL's multidimensional NetCDF driver (so no netcdf4 / h5netcdf engine is required). Reading back with NetCDF.read_file is pure GDAL.
• It is eager: it materialises the full (time, band, y, x) cube in memory. Above ~2 GB it warns; for large cubes prefer DatasetCollection.to_zarr, which streams the data chunk-by-chunk.
• The output is CF-1.8 with the geobox attached (crs_wkt / GeoTransform); the time axis is CF-encoded as nanoseconds since 1970-01-01 and both the calendar axis and the no-data value round-trip through NetCDF.read_file.
• See also: Converting between GeoTIFF and NetCDF and the Zarr examples for very large cubes.