Exporting an ImageCollection to NetCDF

Introduction

XEE is a XArray extension that allows reading data from Google Earth Engine (GEE). This extension is useful for extracting subsets of data hosted in GEE. In this tutorial, we will learn how to extract a large ImageCollection and save it as a local NetCDF file.

Overview of the Task

We will take the ERA5-Land Daily Aggregated collection and export the subset for the chosen country as a CF-compatible format NetCDF file.

Input Layers:

• ne_10m_admin_0_countries_ind.zip: A shapefile of country boundaries

Output:

• era5_subset.nc: A NetCDF subset of the ERA5-Land Daily Aggregated data.

Data Credit:

• ERA5-Land monthly averaged data from 1981 to present. Copernicus Climate Change Service (C3S) Climate Data Store (CDS).

• Made with Natural Earth. Free vector and raster map data @ naturalearthdata.com.

Setup and Data Download

The following blocks of code will install the required packages and download the datasets to your Colab environment.

%%capture
if 'google.colab' in str(get_ipython()):
    !pip install xee rioxarray netCDF4

import ee
import geopandas as gpd
import numpy as np
import os
import rioxarray as rxr
import xarray as xr

data_folder = 'data'
output_folder = 'output'

if not os.path.exists(data_folder):
    os.mkdir(data_folder)
if not os.path.exists(output_folder):
    os.mkdir(output_folder)

def download(url):
    filename = os.path.join(data_folder, os.path.basename(url))
    if not os.path.exists(filename):
        from urllib.request import urlretrieve
        local, _ = urlretrieve(url, filename)
        print('Downloaded ' + local)

shapefile = 'ne_10m_admin_0_countries_ind.zip'

data_url = 'https://github.com/spatialthoughts/geopython-tutorials/releases/download/data/'
download('{}/{}'.format(data_url,shapefile))

Downloaded data/ne_10m_admin_0_countries_ind.zip

Initialize EE with the High-Volume EndPoint recommended to be used with XEE. Replace the value of the cloud_project variable with your own project id that is linked with GEE.

cloud_project = 'spatialthoughts' # replace with your project id

try:
    ee.Initialize(
        project=cloud_project,
        opt_url='https://earthengine-highvolume.googleapis.com')
except:
    ee.Authenticate()
    ee.Initialize(
        project=cloud_project,
        opt_url='https://earthengine-highvolume.googleapis.com')

Data Preparation

We read the Natural Earth administrative regions shapefile and select a country.

shapefile_path = os.path.join(data_folder, shapefile)
boundaries_gdf = gpd.read_file(shapefile_path)

Select the country of your choice using the 3-letter ISO code for your country. Here we use GHA for Ghana.

country = boundaries_gdf[boundaries_gdf['ADM0_A3'] == 'GHA']
geometry = country.geometry.union_all()

Configure the time period and variables.

start_year = 2023
end_year = 2024
variable = 'temperature_2m'

Open the subset using XEE. Define the ImageCollection and apply filters using the Earth Engine Python API syntax.

start_date = ee.Date.fromYMD(start_year, 1, 1)
end_date = ee.Date.fromYMD(end_year + 1, 1, 1)
bbox = ee.Algorithms.GeometryConstructors.BBox(*geometry.bounds)

era5 = ee.ImageCollection('ECMWF/ERA5_LAND/DAILY_AGGR')
filtered = era5 \
  .filter(ee.Filter.date(start_date, end_date)) \
  .filter(ee.Filter.bounds(bbox)) \
  .select(variable)

XArray can open the filtered collection directly from GEE if the XEE extension is installed.

ds = xr.open_dataset(
    filtered,
    engine='ee',
    geometry=bbox,
    projection=filtered.first().select(0).projection(),
    fast_time_slicing=True,
    chunks={'index': 24},
)
ds

<xarray.Dataset> Size: 8MB
Dimensions:         (time: 731, lon: 45, lat: 64)
Coordinates:
  * time            (time) datetime64[ns] 6kB 2023-01-01 ... 2024-12-31
  * lon             (lon) float64 360B -3.213 -3.113 -3.013 ... 1.087 1.187
  * lat             (lat) float64 512B 11.11 11.01 10.91 ... 5.013 4.913 4.813
Data variables:
    temperature_2m  (time, lon, lat) float32 8MB dask.array<chunksize=(48, 45, 64), meta=np.ndarray>
Attributes:
    crs:      EPSG:4326

xarray.Dataset

• Dimensions:
  • time: 731
  • lon: 45
  • lat: 64
• Coordinates: (3)
  • time
    (time)
    datetime64[ns]
    2023-01-01 ... 2024-12-31

    array(['2023-01-01T00:00:00.000000000', '2023-01-02T00:00:00.000000000',
           '2023-01-03T00:00:00.000000000', ..., '2024-12-29T00:00:00.000000000',
           '2024-12-30T00:00:00.000000000', '2024-12-31T00:00:00.000000000'],
          dtype='datetime64[ns]')

  • lon
    (lon)
    float64
    -3.213 -3.113 ... 1.087 1.187

    array([-3.212509, -3.112509, -3.012509, -2.912509, -2.812509, -2.712509,
           -2.612509, -2.512509, -2.412509, -2.312509, -2.212509, -2.112509,
           -2.012509, -1.912509, -1.812509, -1.712509, -1.612509, -1.512509,
           -1.412509, -1.312509, -1.212509, -1.112509, -1.012509, -0.912509,
           -0.812509, -0.712509, -0.612509, -0.512509, -0.412509, -0.312509,
           -0.212509, -0.112509, -0.012509,  0.087491,  0.187491,  0.287491,
            0.387491,  0.487491,  0.587491,  0.687491,  0.787491,  0.887491,
            0.987491,  1.087491,  1.187491])

  • lat
    (lat)
    float64
    11.11 11.01 10.91 ... 4.913 4.813

    array([11.112937, 11.012937, 10.912937, 10.812937, 10.712937, 10.612937,
           10.512937, 10.412937, 10.312937, 10.212937, 10.112937, 10.012937,
            9.912937,  9.812937,  9.712937,  9.612937,  9.512937,  9.412937,
            9.312937,  9.212937,  9.112937,  9.012937,  8.912937,  8.812937,
            8.712937,  8.612937,  8.512937,  8.412937,  8.312937,  8.212937,
            8.112937,  8.012937,  7.912937,  7.812937,  7.712937,  7.612937,
            7.512937,  7.412937,  7.312937,  7.212937,  7.112937,  7.012937,
            6.912937,  6.812937,  6.712937,  6.612937,  6.512937,  6.412937,
            6.312937,  6.212937,  6.112937,  6.012937,  5.912937,  5.812937,
            5.712937,  5.612937,  5.512937,  5.412937,  5.312937,  5.212937,
            5.112937,  5.012937,  4.912937,  4.812937])

• Data variables: (1)
  • temperature_2m
    (time, lon, lat)
    float32
    dask.array<chunksize=(48, 45, 64), meta=np.ndarray>

    id :

    temperature_2m

    data_type :

    {'type': 'PixelType', 'precision': 'double'}

    dimensions :

    [3601, 1801]

    crs :

    EPSG:4326

    crs_transform :

    [0.1, 0, -180.05, 0, -0.1, 90.05]

    Array Chunk Bytes 8.03 MiB 540.00 kiB Shape (731, 45, 64) (48, 45, 64) Dask graph 16 chunks in 2 graph layers Data type float32 numpy.ndarray

• Indexes: (3)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2023-01-01', '2023-01-02', '2023-01-03', '2023-01-04',
                   '2023-01-05', '2023-01-06', '2023-01-07', '2023-01-08',
                   '2023-01-09', '2023-01-10',
                   ...
                   '2024-12-22', '2024-12-23', '2024-12-24', '2024-12-25',
                   '2024-12-26', '2024-12-27', '2024-12-28', '2024-12-29',
                   '2024-12-30', '2024-12-31'],
                  dtype='datetime64[ns]', name='time', length=731, freq=None))

  • lon
    PandasIndex

    PandasIndex(Index([  -3.2125093179999853,    -3.112509317999985,    -3.012509317999985,
              -2.912509317999985,    -2.812509317999985,   -2.7125093179999853,
              -2.612509317999985,    -2.512509317999985,    -2.412509317999985,
              -2.312509317999985,   -2.2125093179999853,   -2.1125093179999848,
              -2.012509317999985,    -1.912509317999985,    -1.812509317999985,
              -1.712509317999985,    -1.612509317999985,   -1.5125093179999851,
              -1.412509317999985,    -1.312509317999985,   -1.2125093179999848,
             -1.1125093179999852,   -1.0125093179999851,    -0.912509317999985,
             -0.8125093179999849,   -0.7125093179999848,   -0.6125093179999848,
             -0.5125093179999851,    -0.412509317999985,  -0.31250931799998494,
            -0.21250931799998485,  -0.11250931799998476, -0.012509317999985115,
             0.08749068200001497,   0.18749068200001506,   0.28749068200001515,
             0.38749068200001524,    0.4874906820000149,     0.587490682000015,
              0.6874906820000151,    0.7874906820000147,    0.8874906820000152,
              0.9874906820000149,    1.0874906820000154,     1.187490682000015],
          dtype='float64', name='lon'))

  • lat
    PandasIndex

    PandasIndex(Index([11.112937317000043, 11.012937317000043, 10.912937317000043,
           10.812937317000044, 10.712937317000044, 10.612937317000043,
           10.512937317000043, 10.412937317000043, 10.312937317000044,
           10.212937317000044, 10.112937317000043, 10.012937317000043,
            9.912937317000043,  9.812937317000044,  9.712937317000044,
            9.612937317000043,  9.512937317000043,  9.412937317000043,
            9.312937317000044,  9.212937317000044,  9.112937317000043,
            9.012937317000043,  8.912937317000043,  8.812937317000044,
            8.712937317000044,  8.612937317000043,  8.512937317000043,
            8.412937317000043,  8.312937317000044,  8.212937317000044,
            8.112937317000043,  8.012937317000043,  7.912937317000043,
            7.812937317000044,  7.712937317000043,  7.612937317000043,
            7.512937317000043,  7.412937317000043,  7.312937317000044,
            7.212937317000043,  7.112937317000044,  7.012937317000043,
            6.912937317000043,  6.812937317000043,  6.712937317000043,
            6.612937317000044,  6.512937317000043,  6.412937317000043,
            6.312937317000043,  6.212937317000043,  6.112937317000043,
            6.012937317000043,  5.912937317000043,  5.812937317000043,
            5.712937317000043,  5.612937317000043,  5.512937317000043,
            5.412937317000043,  5.312937317000043,  5.212937317000043,
            5.112937317000043,  5.012937317000043,  4.912937317000043,
            4.812937317000043],
          dtype='float64', name='lat'))

• Attributes: (1)

  crs :

  EPSG:4326

XEE queries the data for the bounding box of the provided geometry. So we still need to clip the pixels outside the geometry. We rename the lat and lon dimensions to y and x to make it compatible with rioxarray and clip it.

clipped_ds = ds \
  .rename({'lon': 'x', 'lat':'y'}) \
  .rio.clip(country.geometry.values)

Export

Prepare the data for export. Currently XEE produces dataset with dimension ordering that is different than what is expected by many other Python/R packages. We transpose the data to fix this as recommended by XEE developers.

export_ds = clipped_ds.transpose()

# Plot a single timestep to verify the data
export_ds[variable].isel(time=-1).plot.imshow()

<matplotlib.image.AxesImage at 0x7e082e560e50>

Save the results as a NetCDF file.

%%time
output_file = f'era5_subset.nc'
output_path = os.path.join(output_folder, output_file)

# Enable compression
encoding = {variable: {'zlib': True}}
export_ds.to_netcdf(output_path, encoding=encoding)

CPU times: user 567 ms, sys: 105 ms, total: 672 ms
Wall time: 6.76 s

---

If you want to give feedback or share your experience with this tutorial, please comment below. (requires GitHub account)