Transform NorESM data to DGGS

Context

Purpose

The goal is to experiment with how to transform NorESM data into DGGS (healpix grid).

Warning

This FAIR2Adapt case study will likely require different types of interpolation, as bilinear interpolation is not suitable for all variables. For instance, conservative interpolation is currently being developed.

Description

In this notebook, we will:

• Open online NorESM data

• Transform one variable (TREFHT) into DGGS Healpix

• Save the transformed data in Zarr

Contributions

Notebook

• Even Moa Myklebust, Simula Research Laboratory (Norway) (author), @evenmm

• Anne Fouilloux, Simula Research Laboratory (Norway) (reviewer), @annefou

• Ola Formo Kihle, Independent Consultant / UW Contractor (Norway) (reviewer), @ofk123

Bibliography and other interesting resources

• The Norwegian Earth System Model (NorESM)

• NorESM dataset used in this notebook, 10.11582/2024.00093

Import Libraries

# Install xarray-healpy and dggs libraries for regridding
#!pip install git+https://github.com/IAOCEA/xarray-healpy.git git+https://github.com/xarray-contrib/xdggs.git cmcrameri

import os
import numpy as np
import s3fs                                                # For accessing file storage at Amazon S3
import xdggs                                               # Discrete global grid systems in x-array
import xarray as xr                                        # N-dimensional arrays with dimension, coordinate and attribute labels
from xarray_healpy import HealpyGridInfo, HealpyRegridder
xr.set_options(display_expand_data=False, display_expand_attrs=False, keep_attrs=True)

<xarray.core.options.set_options at 0x7813af5349b0>

Open file with Xarray

• The file is on the cloud (accessible to everyone)

• The first time, opening the atmosphere ds typically takes a few minutes

• Therefore, we save a local version with selected subset of variables for faster loading the nezxt times

# Define file paths
subset_dataset_path = "small_N1850_f19_tn14_20190722.cam.h0.1801-01.zarr"
endpoint_url = 'https://pangeo-eosc-minioapi.vm.fedcloud.eu/'
filename_atmosphere = "s3://afouilloux-fair2adapt/10.11582_2024.00093/N1850_f19_tn14_20190722/atm/hist/N1850_f19_tn14_20190722.cam.h0.1801-01.nc"

# Check if the small dataset exists
if os.path.exists(subset_dataset_path):
    print(f"Loading local subset dataset from {subset_dataset_path}")
    ds = xr.open_zarr(subset_dataset_path)
else:
    # Extract file from S3
    print(f"Local subset of variables not found. Loading full dataset from {endpoint_url}")
    client_kwargs={'endpoint_url': endpoint_url}
    s3 = s3fs.S3FileSystem(anon=True, client_kwargs=client_kwargs)
    
    # Opening the atmosphere ds typically takes a few minutes
    atmosphere_ds = xr.open_dataset(s3.open(filename_atmosphere))

    # Save a subset for faster access on reload
    print("Saving subset of variables to subset_dataset")
    atmosphere_ds[["SST", "TREFHT"]].to_zarr(subset_dataset_path)
    ds = xr.open_zarr(subset_dataset_path)

# Display the subset dataset
ds

Loading local subset dataset from small_N1850_f19_tn14_20190722.cam.h0.1801-01.zarr

<xarray.Dataset> Size: 113kB
Dimensions:  (time: 1, lat: 96, lon: 144)
Coordinates:
  * lat      (lat) float64 768B -90.0 -88.11 -86.21 -84.32 ... 86.21 88.11 90.0
  * lon      (lon) float64 1kB 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
  * time     (time) object 8B 1801-02-01 00:00:00
Data variables:
    SST      (time, lat, lon) float32 55kB dask.array<chunksize=(1, 96, 144), meta=np.ndarray>
    TREFHT   (time, lat, lon) float32 55kB dask.array<chunksize=(1, 96, 144), meta=np.ndarray>
Attributes: (9)

xarray.Dataset

• Dimensions:
  • time: 1
  • lat: 96
  • lon: 144
• Coordinates: (3)
  • lat
    (lat)
    float64
    -90.0 -88.11 -86.21 ... 88.11 90.0

    long_name :

    latitude

    units :

    degrees_north

    array([-90.      , -88.105263, -86.210526, -84.315789, -82.421053, -80.526316,
           -78.631579, -76.736842, -74.842105, -72.947368, -71.052632, -69.157895,
           -67.263158, -65.368421, -63.473684, -61.578947, -59.684211, -57.789474,
           -55.894737, -54.      , -52.105263, -50.210526, -48.315789, -46.421053,
           -44.526316, -42.631579, -40.736842, -38.842105, -36.947368, -35.052632,
           -33.157895, -31.263158, -29.368421, -27.473684, -25.578947, -23.684211,
           -21.789474, -19.894737, -18.      , -16.105263, -14.210526, -12.315789,
           -10.421053,  -8.526316,  -6.631579,  -4.736842,  -2.842105,  -0.947368,
             0.947368,   2.842105,   4.736842,   6.631579,   8.526316,  10.421053,
            12.315789,  14.210526,  16.105263,  18.      ,  19.894737,  21.789474,
            23.684211,  25.578947,  27.473684,  29.368421,  31.263158,  33.157895,
            35.052632,  36.947368,  38.842105,  40.736842,  42.631579,  44.526316,
            46.421053,  48.315789,  50.210526,  52.105263,  54.      ,  55.894737,
            57.789474,  59.684211,  61.578947,  63.473684,  65.368421,  67.263158,
            69.157895,  71.052632,  72.947368,  74.842105,  76.736842,  78.631579,
            80.526316,  82.421053,  84.315789,  86.210526,  88.105263,  90.      ])

  • lon
    (lon)
    float64
    0.0 2.5 5.0 ... 352.5 355.0 357.5

    long_name :

    longitude

    units :

    degrees_east

    array([  0. ,   2.5,   5. ,   7.5,  10. ,  12.5,  15. ,  17.5,  20. ,  22.5,
            25. ,  27.5,  30. ,  32.5,  35. ,  37.5,  40. ,  42.5,  45. ,  47.5,
            50. ,  52.5,  55. ,  57.5,  60. ,  62.5,  65. ,  67.5,  70. ,  72.5,
            75. ,  77.5,  80. ,  82.5,  85. ,  87.5,  90. ,  92.5,  95. ,  97.5,
           100. , 102.5, 105. , 107.5, 110. , 112.5, 115. , 117.5, 120. , 122.5,
           125. , 127.5, 130. , 132.5, 135. , 137.5, 140. , 142.5, 145. , 147.5,
           150. , 152.5, 155. , 157.5, 160. , 162.5, 165. , 167.5, 170. , 172.5,
           175. , 177.5, 180. , 182.5, 185. , 187.5, 190. , 192.5, 195. , 197.5,
           200. , 202.5, 205. , 207.5, 210. , 212.5, 215. , 217.5, 220. , 222.5,
           225. , 227.5, 230. , 232.5, 235. , 237.5, 240. , 242.5, 245. , 247.5,
           250. , 252.5, 255. , 257.5, 260. , 262.5, 265. , 267.5, 270. , 272.5,
           275. , 277.5, 280. , 282.5, 285. , 287.5, 290. , 292.5, 295. , 297.5,
           300. , 302.5, 305. , 307.5, 310. , 312.5, 315. , 317.5, 320. , 322.5,
           325. , 327.5, 330. , 332.5, 335. , 337.5, 340. , 342.5, 345. , 347.5,
           350. , 352.5, 355. , 357.5])

  • time
    (time)
    object
    1801-02-01 00:00:00

    bounds :

    time_bnds

    long_name :

    time

    array([cftime.DatetimeNoLeap(1801, 2, 1, 0, 0, 0, 0, has_year_zero=True)],
          dtype=object)

• Data variables: (2)
  • SST
    (time, lat, lon)
    float32
    dask.array<chunksize=(1, 96, 144), meta=np.ndarray>

    cell_methods :

    time: mean

    long_name :

    sea surface temperature

    units :

    K

    Array Chunk Bytes 54.00 kiB 54.00 kiB Shape (1, 96, 144) (1, 96, 144) Dask graph 1 chunks in 2 graph layers Data type float32 numpy.ndarray

  • TREFHT
    (time, lat, lon)
    float32
    dask.array<chunksize=(1, 96, 144), meta=np.ndarray>

    cell_methods :

    time: mean

    long_name :

    Reference height temperature

    units :

    K

    Array Chunk Bytes 54.00 kiB 54.00 kiB Shape (1, 96, 144) (1, 96, 144) Dask graph 1 chunks in 2 graph layers Data type float32 numpy.ndarray

• Indexes: (3)
  • lat
    PandasIndex

    PandasIndex(Index([              -90.0,  -88.10526315789474,  -86.21052631578948,
             -84.3157894736842,  -82.42105263157895,  -80.52631578947368,
            -78.63157894736842,  -76.73684210526316,  -74.84210526315789,
            -72.94736842105263,  -71.05263157894737,  -69.15789473684211,
            -67.26315789473685,  -65.36842105263158, -63.473684210526315,
           -61.578947368421055, -59.684210526315795,  -57.78947368421053,
            -55.89473684210527,               -54.0,  -52.10526315789474,
            -50.21052631578947,  -48.31578947368421,  -46.42105263157895,
           -44.526315789473685, -42.631578947368425,  -40.73684210526316,
             -38.8421052631579,  -36.94736842105264,  -35.05263157894737,
            -33.15789473684211, -31.263157894736842, -29.368421052631582,
           -27.473684210526322, -25.578947368421055, -23.684210526315795,
           -21.789473684210535,  -19.89473684210526,               -18.0,
            -16.10526315789474,  -14.21052631578948,  -12.31578947368422,
           -10.421052631578945,  -8.526315789473685,  -6.631578947368425,
            -4.736842105263165, -2.8421052631579045, -0.9473684210526301,
            0.9473684210526301,  2.8421052631578902,    4.73684210526315,
             6.631578947368411,   8.526315789473685,  10.421052631578945,
            12.315789473684205,  14.210526315789465,  16.105263157894726,
                          18.0,   19.89473684210526,   21.78947368421052,
             23.68421052631578,   25.57894736842104,  27.473684210526315,
            29.368421052631575,  31.263157894736835,  33.157894736842096,
            35.052631578947356,   36.94736842105263,   38.84210526315789,
             40.73684210526315,   42.63157894736841,   44.52631578947367,
             46.42105263157893,   48.31578947368419,   50.21052631578948,
             52.10526315789474,                54.0,   55.89473684210526,
             57.78947368421052,   59.68421052631578,   61.57894736842104,
              63.4736842105263,   65.36842105263156,   67.26315789473682,
             69.15789473684211,   71.05263157894737,   72.94736842105263,
             74.84210526315789,   76.73684210526315,   78.63157894736841,
             80.52631578947367,   82.42105263157893,   84.31578947368419,
             86.21052631578945,   88.10526315789474,                90.0],
          dtype='float64', name='lat'))

  • lon
    PandasIndex

    PandasIndex(Index([  0.0,   2.5,   5.0,   7.5,  10.0,  12.5,  15.0,  17.5,  20.0,  22.5,
           ...
           335.0, 337.5, 340.0, 342.5, 345.0, 347.5, 350.0, 352.5, 355.0, 357.5],
          dtype='float64', name='lon', length=144))

  • time
    PandasIndex

    PandasIndex(CFTimeIndex([1801-02-01 00:00:00],
                dtype='object', length=1, calendar='noleap', freq=None))

• Attributes: (9)

  Conventions :

  CF-1.0

  case :

  N1850_f19_tn14_20190722

  host :

  initial_file :

  /cluster/shared/noresm/inputdata/atm/cam/inic/fv/cami-mam3_0000-01-01_1.9x2.5_L32_c150407.nc

  logname :

  olivie

  model_doi_url :

  https://doi.org/10.5065/D67H1H0V

  source :

  CAM

  time_period_freq :

  month_1

  topography_file :

  /cluster/shared/noresm/inputdata/atm/cam/topo/fv_1.9x2.5_nc3000_Nsw084_Nrs016_Co120_Fi001_ZR_061116.nc

Quick visualization of reference height temperature (TREFHT)

# Plot reference height temperature (TREFHT)
ds.TREFHT.isel(time=0).plot()

<matplotlib.collections.QuadMesh at 0x78133dbe8260>

Standardize coordinate names and center longitude

ds = ds.rename_dims({"lat": "latitude", "lon": "longitude"})
ds = ds.rename_vars({"lat": "latitude", "lon": "longitude"})
ds.latitude.attrs["standard_name"] = 'latitude'
ds.longitude.attrs["standard_name"] = 'longitude'

Center around 0 longitude:

ds.coords['longitude'] = (ds.coords['longitude'] + 180) % 360 - 180
ds = ds.sortby(ds.longitude)

Regridding

Define the target grid (Healpix DGGS)

def center_longitude(ds, center):
    if center == 0:
        centered = (ds.longitude + 180) % 360 - 180
    else:
        centered = (ds.longitude - 180) % 360 + 180
    return ds.assign_coords(longitude=centered)

nside = 16                               # Each side of the original 12 faces in Healpix is divided into nside parts 
healpy_grid_level = int(np.log2(nside))  # Healpix level
number_of_cells = 12*nside**2            # The resulting total number of cells

print("nside:", nside)
print("Level:", healpy_grid_level)
print("Number of cells:", number_of_cells)

nside: 16
Level: 4
Number of cells: 3072

# Define the target Healpix grid
grid = HealpyGridInfo(level=healpy_grid_level)
target_grid = grid.target_grid(ds).pipe(center_longitude, 0).drop_attrs(deep=False).dggs.decode({"grid_name": "healpix", "nside": nside, "indexing_scheme": "nested"})

%%time
# Compute interpolation weights for regridding diff data
regridder = HealpyRegridder(
    ds[["longitude", "latitude"]].compute(),
    target_grid,
    method="bilinear",
    interpolation_kwargs={},
)
print(regridder)

HealpyRegridder(method='bilinear', interpolation_kwargs={})
CPU times: user 3.59 s, sys: 228 ms, total: 3.82 s
Wall time: 3.83 s

%%time
regridded = regridder.regrid_ds(ds).pipe(xdggs.decode)

CPU times: user 1.38 s, sys: 616 ms, total: 1.99 s
Wall time: 2.01 s

%%time
ds_regridded = regridded.TREFHT.compute().squeeze()

CPU times: user 186 ms, sys: 12.1 ms, total: 198 ms
Wall time: 196 ms

ds_regridded.dggs.explore()

Save the transformed data to zarr

url = f"./NorESM-TREFHT-healpix-lvl-{healpy_grid_level}.zarr"
print("url:", url)
ds_regridded.to_zarr(url,mode='w')

url: ./NorESM-TREFHT-healpix-lvl-4.zarr

<xarray.backends.zarr.ZarrStore at 0x78132be8eec0>