Vertical Regridding

Authors: Jason Boutte and Jill Zhang

Date: 9/26/23

Related APIs:

• xarray.Dataset.regridder.vertical

The data used in this example can be found through the Earth System Grid Federation (ESGF) search portal. We are using xarray’s OPeNDAP support to read a netCDF4 dataset file directly from its source. The data is not loaded over the network until we perform operations on it (e.g., temperature unit adjustment). More information on the xarray’s OPeNDAP support can be found here.

We’ll cover vertical regridding using xgcm. Two examples are outlined here to apply vertical regridding/remapping using ocean variables and atmosphere variables, respectively.

Notebook Setup

Create an Anaconda environment for this notebook using the command below, then select the kernel in Jupyter.

conda create -n xcdat-notebook -c conda-forge python xarray netcdf4 xcdat xesmf matplotlib nc-time-axis jupyter

• xesmf is required for horizontal regridding with xESMF

• matplotlib is an optional dependency required for plotting with xarray

• nc-time-axis is an optional dependency required for matplotlib to plot cftime coordinates

Example 1: Remapping Ocean Variables

%matplotlib inline

import matplotlib.pyplot as plt
import xarray as xr
import xcdat
import numpy as np
# gsw_xarray is a wrapper for GSW-Python:
# the Python implementation of the Gibbs SeaWater (GSW) Oceanographic Toolbox of TEOS-10
import gsw_xarray as gsw

import warnings

warnings.filterwarnings("ignore")

1. Open dataset

# urls for sea water potential temperature (thetao) and salinity (so) from the NCAR model in CMIP6
urls = [
    "http://aims3.llnl.gov/thredds/dodsC/css03_data/CMIP6/CMIP/NCAR/CESM2/historical/r1i1p1f1/Omon/so/gn/v20190308/so_Omon_CESM2_historical_r1i1p1f1_gn_185001-201412.nc",
    "http://aims3.llnl.gov/thredds/dodsC/css03_data/CMIP6/CMIP/NCAR/CESM2/historical/r1i1p1f1/Omon/thetao/gn/v20190308/thetao_Omon_CESM2_historical_r1i1p1f1_gn_185001-201412.nc",
]

ds = xr.merge([xcdat.open_dataset(x, chunks={"time": 4}) for x in urls])

# lev coordinate is in cm and bounds is in m, convert lev to m
with xr.set_options(keep_attrs=True):
    ds.lev.load()
    ds["lev"] = ds.lev / 100
    ds.lev.attrs["units"] = "meters"

ds

2023-09-26 15:45:48,860 [WARNING]: bounds.py(add_missing_bounds:186) >> The nlat coord variable has a 'units' attribute that is not in degrees.
2023-09-26 15:45:48,860 [WARNING]: bounds.py(add_missing_bounds:186) >> The nlat coord variable has a 'units' attribute that is not in degrees.
2023-09-26 15:45:48,966 [WARNING]: bounds.py(add_missing_bounds:186) >> The nlat coord variable has a 'units' attribute that is not in degrees.
2023-09-26 15:45:48,966 [WARNING]: bounds.py(add_missing_bounds:186) >> The nlat coord variable has a 'units' attribute that is not in degrees.

<xarray.Dataset>
Dimensions:    (lev: 60, nlat: 384, nlon: 320, time: 1980, d2: 2, vertices: 4,
                bnds: 2)
Coordinates:
  * lev        (lev) float64 5.0 15.0 25.0 ... 4.875e+03 5.125e+03 5.375e+03
  * nlat       (nlat) int32 1 2 3 4 5 6 7 8 ... 377 378 379 380 381 382 383 384
  * nlon       (nlon) int32 1 2 3 4 5 6 7 8 ... 313 314 315 316 317 318 319 320
  * time       (time) object 1850-01-15 13:00:00.000007 ... 2014-12-15 12:00:00
    lat        (nlat, nlon) float64 dask.array<chunksize=(384, 320), meta=np.ndarray>
    lon        (nlat, nlon) float64 dask.array<chunksize=(384, 320), meta=np.ndarray>
Dimensions without coordinates: d2, vertices, bnds
Data variables:
    time_bnds  (time, d2) object dask.array<chunksize=(4, 2), meta=np.ndarray>
    lat_bnds   (nlat, nlon, vertices) float32 dask.array<chunksize=(384, 320, 4), meta=np.ndarray>
    lon_bnds   (nlat, nlon, vertices) float32 dask.array<chunksize=(384, 320, 4), meta=np.ndarray>
    lev_bnds   (lev, d2) float32 dask.array<chunksize=(60, 2), meta=np.ndarray>
    so         (time, lev, nlat, nlon) float32 dask.array<chunksize=(4, 60, 384, 320), meta=np.ndarray>
    nlon_bnds  (nlon, bnds) float64 0.5 1.5 1.5 2.5 ... 318.5 319.5 319.5 320.5
    thetao     (time, lev, nlat, nlon) float32 dask.array<chunksize=(4, 60, 384, 320), meta=np.ndarray>
Attributes: (12/46)
    Conventions:                     CF-1.7 CMIP-6.2
    activity_id:                     CMIP
    case_id:                         15
    cesm_casename:                   b.e21.BHIST.f09_g17.CMIP6-historical.001
    contact:                         cesm_cmip6@ucar.edu
    creation_date:                   2019-01-16T23:15:40Z
    ...                              ...
    sub_experiment_id:               none
    branch_time_in_parent:           219000.0
    branch_time_in_child:            674885.0
    branch_method:                   standard
    further_info_url:                https://furtherinfo.es-doc.org/CMIP6.NCA...
    DODS_EXTRA.Unlimited_Dimension:  time

xarray.Dataset

• Dimensions:
  • lev: 60
  • nlat: 384
  • nlon: 320
  • time: 1980
  • d2: 2
  • vertices: 4
  • bnds: 2
• Coordinates: (6)
  • lev
    (lev)
    float64
    5.0 15.0 ... 5.125e+03 5.375e+03

    axis :

    Z

    bounds :

    lev_bnds

    positive :

    down

    standard_name :

    olevel

    title :

    ocean model level

    type :

    double

    units :

    meters

    _ChunkSizes :

    60

    array([5.000000e+00, 1.500000e+01, 2.500000e+01, 3.500000e+01, 4.500000e+01,
           5.500000e+01, 6.500000e+01, 7.500000e+01, 8.500000e+01, 9.500000e+01,
           1.050000e+02, 1.150000e+02, 1.250000e+02, 1.350000e+02, 1.450000e+02,
           1.550000e+02, 1.650984e+02, 1.754790e+02, 1.862913e+02, 1.976603e+02,
           2.097114e+02, 2.225783e+02, 2.364088e+02, 2.513702e+02, 2.676542e+02,
           2.854837e+02, 3.051192e+02, 3.268680e+02, 3.510935e+02, 3.782276e+02,
           4.087846e+02, 4.433777e+02, 4.827367e+02, 5.277280e+02, 5.793729e+02,
           6.388626e+02, 7.075633e+02, 7.870025e+02, 8.788252e+02, 9.847059e+02,
           1.106204e+03, 1.244567e+03, 1.400497e+03, 1.573946e+03, 1.764003e+03,
           1.968944e+03, 2.186457e+03, 2.413972e+03, 2.649001e+03, 2.889385e+03,
           3.133405e+03, 3.379793e+03, 3.627670e+03, 3.876452e+03, 4.125768e+03,
           4.375392e+03, 4.625190e+03, 4.875083e+03, 5.125028e+03, 5.375000e+03])

  • nlat
    (nlat)
    int32
    1 2 3 4 5 6 ... 380 381 382 383 384

    long_name :

    cell index along second dimension

    units :

    1

    _ChunkSizes :

    384

    array([  1,   2,   3, ..., 382, 383, 384], dtype=int32)

  • nlon
    (nlon)
    int32
    1 2 3 4 5 6 ... 316 317 318 319 320

    long_name :

    cell index along first dimension

    units :

    1

    _ChunkSizes :

    320

    bounds :

    nlon_bnds

    array([  1,   2,   3, ..., 318, 319, 320], dtype=int32)

  • time
    (time)
    object
    1850-01-15 13:00:00.000007 ... 2...

    axis :

    T

    bounds :

    time_bnds

    standard_name :

    time

    title :

    time

    type :

    double

    _ChunkSizes :

    512

    array([cftime.DatetimeNoLeap(1850, 1, 15, 13, 0, 0, 7, has_year_zero=True),
           cftime.DatetimeNoLeap(1850, 2, 14, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1850, 3, 15, 12, 0, 0, 0, has_year_zero=True),
           ...,
           cftime.DatetimeNoLeap(2014, 10, 15, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2014, 11, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(2014, 12, 15, 12, 0, 0, 0, has_year_zero=True)],
          dtype=object)

  • lat
    (nlat, nlon)
    float64
    dask.array<chunksize=(384, 320), meta=np.ndarray>

    axis :

    Y

    bounds :

    lat_bnds

    standard_name :

    latitude

    title :

    Latitude

    type :

    double

    units :

    degrees_north

    valid_max :

    90.0

    valid_min :

    -90.0

    _ChunkSizes :

    [384 320]

    Array Chunk Bytes 0.94 MiB 0.94 MiB Shape (384, 320) (384, 320) Dask graph 1 chunks in 5 graph layers Data type float64 numpy.ndarray

  • lon
    (nlat, nlon)
    float64
    dask.array<chunksize=(384, 320), meta=np.ndarray>

    axis :

    X

    bounds :

    lon_bnds

    standard_name :

    longitude

    title :

    Longitude

    type :

    double

    units :

    degrees_east

    valid_max :

    360.0

    valid_min :

    0.0

    _ChunkSizes :

    [384 320]

    Array Chunk Bytes 0.94 MiB 0.94 MiB Shape (384, 320) (384, 320) Dask graph 1 chunks in 5 graph layers Data type float64 numpy.ndarray

• Data variables: (7)
  • time_bnds
    (time, d2)
    object
    dask.array<chunksize=(4, 2), meta=np.ndarray>

    _ChunkSizes :

    [1 2]

    Array Chunk Bytes 30.94 kiB 64 B Shape (1980, 2) (4, 2) Dask graph 495 chunks in 8 graph layers Data type object numpy.ndarray

  • lat_bnds
    (nlat, nlon, vertices)
    float32
    dask.array<chunksize=(384, 320, 4), meta=np.ndarray>

    units :

    degrees_north

    _ChunkSizes :

    [384 320 4]

    Array Chunk Bytes 1.88 MiB 1.88 MiB Shape (384, 320, 4) (384, 320, 4) Dask graph 1 chunks in 5 graph layers Data type float32 numpy.ndarray

  • lon_bnds
    (nlat, nlon, vertices)
    float32
    dask.array<chunksize=(384, 320, 4), meta=np.ndarray>

    units :

    degrees_east

    _ChunkSizes :

    [384 320 4]

    Array Chunk Bytes 1.88 MiB 1.88 MiB Shape (384, 320, 4) (384, 320, 4) Dask graph 1 chunks in 5 graph layers Data type float32 numpy.ndarray

  • lev_bnds
    (lev, d2)
    float32
    dask.array<chunksize=(60, 2), meta=np.ndarray>

    units :

    m

    _ChunkSizes :

    [60 2]

    Array Chunk Bytes 480 B 480 B Shape (60, 2) (60, 2) Dask graph 1 chunks in 5 graph layers Data type float32 numpy.ndarray

  • so
    (time, lev, nlat, nlon)
    float32
    dask.array<chunksize=(4, 60, 384, 320), meta=np.ndarray>

    cell_measures :

    area: areacello volume: volcello

    cell_methods :

    area: mean where sea time: mean

    comment :

    Sea water salinity is the salt content of sea water, often on the Practical Salinity Scale of 1978. However, the unqualified term 'salinity' is generic and does not necessarily imply any particular method of calculation. The units of salinity are dimensionless and the units attribute should normally be given as 1e-3 or 0.001 i.e. parts per thousand.

    description :

    Sea water salinity is the salt content of sea water, often on the Practical Salinity Scale of 1978. However, the unqualified term 'salinity' is generic and does not necessarily imply any particular method of calculation. The units of salinity are dimensionless and the units attribute should normally be given as 1e-3 or 0.001 i.e. parts per thousand.

    frequency :

    mon

    id :

    so

    long_name :

    Sea Water Salinity

    mipTable :

    Omon

    out_name :

    so

    prov :

    Omon ((isd.003))

    realm :

    ocean

    standard_name :

    sea_water_salinity

    time :

    time

    time_label :

    time-mean

    time_title :

    Temporal mean

    title :

    Sea Water Salinity

    type :

    real

    units :

    0.001

    variable_id :

    so

    _ChunkSizes :

    [ 1 30 192 160]

    Array Chunk Bytes 54.38 GiB 112.50 MiB Shape (1980, 60, 384, 320) (4, 60, 384, 320) Dask graph 495 chunks in 2 graph layers Data type float32 numpy.ndarray

  • nlon_bnds
    (nlon, bnds)
    float64
    0.5 1.5 1.5 ... 319.5 319.5 320.5

    xcdat_bounds :

    True

    array([[  0.5,   1.5],
           [  1.5,   2.5],
           [  2.5,   3.5],
           [  3.5,   4.5],
           [  4.5,   5.5],
           [  5.5,   6.5],
           [  6.5,   7.5],
           [  7.5,   8.5],
           [  8.5,   9.5],
           [  9.5,  10.5],
           [ 10.5,  11.5],
           [ 11.5,  12.5],
           [ 12.5,  13.5],
           [ 13.5,  14.5],
           [ 14.5,  15.5],
           [ 15.5,  16.5],
           [ 16.5,  17.5],
           [ 17.5,  18.5],
           [ 18.5,  19.5],
           [ 19.5,  20.5],
    ...
           [300.5, 301.5],
           [301.5, 302.5],
           [302.5, 303.5],
           [303.5, 304.5],
           [304.5, 305.5],
           [305.5, 306.5],
           [306.5, 307.5],
           [307.5, 308.5],
           [308.5, 309.5],
           [309.5, 310.5],
           [310.5, 311.5],
           [311.5, 312.5],
           [312.5, 313.5],
           [313.5, 314.5],
           [314.5, 315.5],
           [315.5, 316.5],
           [316.5, 317.5],
           [317.5, 318.5],
           [318.5, 319.5],
           [319.5, 320.5]])

  • thetao
    (time, lev, nlat, nlon)
    float32
    dask.array<chunksize=(4, 60, 384, 320), meta=np.ndarray>

    cell_measures :

    area: areacello volume: volcello

    cell_methods :

    area: mean where sea time: mean

    comment :

    Diagnostic should be contributed even for models using conservative temperature as prognostic field.

    description :

    Diagnostic should be contributed even for models using conservative temperature as prognostic field.

    frequency :

    mon

    id :

    thetao

    long_name :

    Sea Water Potential Temperature

    mipTable :

    Omon

    out_name :

    thetao

    prov :

    Omon ((isd.003))

    realm :

    ocean

    standard_name :

    sea_water_potential_temperature

    time :

    time

    time_label :

    time-mean

    time_title :

    Temporal mean

    title :

    Sea Water Potential Temperature

    type :

    real

    units :

    degC

    variable_id :

    thetao

    _ChunkSizes :

    [ 1 30 192 160]

    Array Chunk Bytes 54.38 GiB 112.50 MiB Shape (1980, 60, 384, 320) (4, 60, 384, 320) Dask graph 495 chunks in 2 graph layers Data type float32 numpy.ndarray

• Indexes: (4)
  • lev
    PandasIndex

    PandasIndex(Index([            5.0,            15.0,            25.0,            35.0,
                      45.0,            55.0,            65.0,            75.0,
                      85.0,            95.0,           105.0,           115.0,
                     125.0,           135.0,           145.0,           155.0,
            165.0983984375, 175.47904296875, 186.29126953125,  197.6602734375,
           209.71138671875,    222.57828125,   236.408828125,    251.37015625,
           267.65419921875, 285.48365234375,    305.11921875, 326.86798828125,
            351.0934765625,  378.2276171875,  408.7846484375,  443.3776953125,
              482.73671875,  527.7280078125,   579.372890625,  638.8626171875,
              707.56328125,        787.0025,   878.825234375,   984.705859375,
             1106.20421875,     1244.566875,    1400.4971875,   1573.94640625,
             1764.00328125,   1968.94421875,    2186.4565625,    2413.9715625,
                2649.00125,    2889.3846875,    3133.4046875,    3379.7934375,
              3627.6703125,     3876.451875,     4125.768125,       4375.3925,
              4625.1903125,    4875.0834375,     5125.028125,          5375.0],
          dtype='float64', name='lev'))

  • nlat
    PandasIndex

    PandasIndex(Index([  1,   2,   3,   4,   5,   6,   7,   8,   9,  10,
           ...
           375, 376, 377, 378, 379, 380, 381, 382, 383, 384],
          dtype='int32', name='nlat', length=384))

  • nlon
    PandasIndex

    PandasIndex(Index([  1,   2,   3,   4,   5,   6,   7,   8,   9,  10,
           ...
           311, 312, 313, 314, 315, 316, 317, 318, 319, 320],
          dtype='int32', name='nlon', length=320))

  • time
    PandasIndex

    PandasIndex(CFTimeIndex([1850-01-15 13:00:00.000007, 1850-02-14 00:00:00, 1850-03-15 12:00:00,
                 1850-04-15 00:00:00, 1850-05-15 12:00:00, 1850-06-15 00:00:00,
                 1850-07-15 12:00:00, 1850-08-15 12:00:00, 1850-09-15 00:00:00,
                 1850-10-15 12:00:00,
                 ...
                 2014-03-15 12:00:00, 2014-04-15 00:00:00, 2014-05-15 12:00:00,
                 2014-06-15 00:00:00, 2014-07-15 12:00:00, 2014-08-15 12:00:00,
                 2014-09-15 00:00:00, 2014-10-15 12:00:00, 2014-11-15 00:00:00,
                 2014-12-15 12:00:00],
                dtype='object', length=1980, calendar='noleap', freq='None'))

• Attributes: (46)

  Conventions :

  CF-1.7 CMIP-6.2

  activity_id :

  CMIP

  case_id :

  15

  cesm_casename :

  b.e21.BHIST.f09_g17.CMIP6-historical.001

  contact :

  cesm_cmip6@ucar.edu

  creation_date :

  2019-01-16T23:15:40Z

  data_specs_version :

  01.00.29

  experiment :

  all-forcing simulation of the recent past

  experiment_id :

  historical

  external_variables :

  areacello volcello

  forcing_index :

  1

  frequency :

  mon

  grid :

  native gx1v7 displaced pole grid (384x320 latxlon)

  grid_label :

  gn

  initialization_index :

  1

  institution :

  National Center for Atmospheric Research, Climate and Global Dynamics Laboratory, 1850 Table Mesa Drive, Boulder, CO 80305, USA

  institution_id :

  NCAR

  license :

  CMIP6 model data produced by <The National Center for Atmospheric Research> is licensed under a Creative Commons Attribution-[]ShareAlike 4.0 International License (https://creativecommons.org/licenses/). Consult https://pcmdi.llnl.gov/CMIP6/TermsOfUse for terms of use governing CMIP6 output, including citation requirements and proper acknowledgment. Further information about this data, including some limitations, can be found via the further_info_url (recorded as a global attribute in this file)[]. The data producers and data providers make no warranty, either express or implied, including, but not limited to, warranties of merchantability and fitness for a particular purpose. All liabilities arising from the supply of the information (including any liability arising in negligence) are excluded to the fullest extent permitted by law.

  mip_era :

  CMIP6

  model_doi_url :

  https://doi.org/10.5065/D67H1H0V

  nominal_resolution :

  100 km

  parent_activity_id :

  CMIP

  parent_experiment_id :

  piControl

  parent_mip_era :

  CMIP6

  parent_source_id :

  CESM2

  parent_time_units :

  days since 0001-01-01 00:00:00

  parent_variant_label :

  r1i1p1f1

  physics_index :

  1

  product :

  model-output

  realization_index :

  1

  realm :

  ocean

  source :

  CESM2 (2017): atmosphere: CAM6 (0.9x1.25 finite volume grid; 288 x 192 longitude/latitude; 32 levels; top level 2.25 mb); ocean: POP2 (320x384 longitude/latitude; 60 levels; top grid cell 0-10 m); sea_ice: CICE5.1 (same grid as ocean); land: CLM5 0.9x1.25 finite volume grid; 288 x 192 longitude/latitude; 32 levels; top level 2.25 mb); aerosol: MAM4 (0.9x1.25 finite volume grid; 288 x 192 longitude/latitude; 32 levels; top level 2.25 mb); atmoschem: MAM4 (0.9x1.25 finite volume grid; 288 x 192 longitude/latitude; 32 levels; top level 2.25 mb); landIce: CISM2.1; ocnBgchem: MARBL (320x384 longitude/latitude; 60 levels; top grid cell 0-10 m)

  source_id :

  CESM2

  source_type :

  AOGCM BGC

  table_id :

  Omon

  tracking_id :

  hdl:21.14100/3fc36959-d797-4d0e-a61a-814cb7e04aa3

  variable_id :

  so

  variant_info :

  CMIP6 20th century experiments (1850-2014) with CAM6, interactive land (CLM5), coupled ocean (POP2) with biogeochemistry (MARBL), interactive sea ice (CICE5.1), and non-evolving land ice (CISM2.1)

  variant_label :

  r1i1p1f1

  sub_experiment :

  none

  sub_experiment_id :

  none

  branch_time_in_parent :

  219000.0

  branch_time_in_child :

  674885.0

  branch_method :

  standard

  further_info_url :

  https://furtherinfo.es-doc.org/CMIP6.NCAR.CESM2.historical.none.r1i1p1f1

  DODS_EXTRA.Unlimited_Dimension :

  time

2. Create the output grid

Related API: xcdat.create_grid()

In this example, we will generate a vertical grid with a linear spaced level coordinate using xcdat.create_grid

Alternatively a grid can be loaded from a file, e.g.

grid_urlpath = "http://aims3.llnl.gov/thredds/dodsC/css03_data/CMIP6/CMIP/NOAA-GFDL/GFDL-CM4/abrupt-4xCO2/r1i1p1f1/day/tas/gr2/v20180701/tas_day_GFDL-CM4_abrupt-4xCO2_r1i1p1f1_gr2_00010101-00201231.nc"

grid_ds = xcdat.open_dataset(grid_urlpath)

output_grid = grid_ds.regridder.grid

output_grid = xcdat.create_grid(
    z=xcdat.create_axis("lev", np.linspace(5, 537, 10))
)

output_grid

<xarray.Dataset>
Dimensions:   (lev: 10, bnds: 2)
Coordinates:
  * lev       (lev) float64 5.0 64.11 123.2 182.3 ... 359.7 418.8 477.9 537.0
Dimensions without coordinates: bnds
Data variables:
    lev_bnds  (lev, bnds) float64 -24.56 34.56 34.56 93.67 ... 507.4 507.4 566.6

xarray.Dataset

• Dimensions:
  • lev: 10
  • bnds: 2
• Coordinates: (1)
  • lev
    (lev)
    float64
    5.0 64.11 123.2 ... 477.9 537.0

    axis :

    Z

    coordinate :

    vertical

    bounds :

    lev_bnds

    array([  5.      ,  64.111111, 123.222222, 182.333333, 241.444444, 300.555556,
           359.666667, 418.777778, 477.888889, 537.      ])

• Data variables: (1)
  • lev_bnds
    (lev, bnds)
    float64
    -24.56 34.56 34.56 ... 507.4 566.6

    xcdat_bounds :

    True

    array([[-24.55555556,  34.55555556],
           [ 34.55555556,  93.66666667],
           [ 93.66666667, 152.77777778],
           [152.77777778, 211.88888889],
           [211.88888889, 271.        ],
           [271.        , 330.11111111],
           [330.11111111, 389.22222222],
           [389.22222222, 448.33333333],
           [448.33333333, 507.44444444],
           [507.44444444, 566.55555556]])

• Indexes: (1)
  • lev
    PandasIndex

    PandasIndex(Index([               5.0,  64.11111111111111, 123.22222222222223,
           182.33333333333334, 241.44444444444446, 300.55555555555554,
            359.6666666666667,  418.7777777777778,  477.8888888888889,
                        537.0],
          dtype='float64', name='lev'))

• Attributes: (0)

3. Regridding using the linear method

Related API: xarray.Dataset.regridder.vertical()

Here we will regrid the input data to the output grid using the xgcm tool and the linear method.

We’ll interpolate salinity onto the new vertical grid.

output = ds.regridder.vertical("so", output_grid, tool="xgcm", method="linear")

output.so.isel(time=0).mean(dim="nlon").plot()
plt.gca().invert_yaxis()

2023-09-26 15:45:50,454 [WARNING]: bounds.py(add_missing_bounds:186) >> The nlat coord variable has a 'units' attribute that is not in degrees.
2023-09-26 15:45:50,454 [WARNING]: bounds.py(add_missing_bounds:186) >> The nlat coord variable has a 'units' attribute that is not in degrees.

4. Regridding from depth to density space

Related API: xarray.Dataset.regridder.vertical()

Here we will regrid the input data to the output grid using the xgcm tool and the linear method.

We’ll remap salinity into density space.

# Apply gsw function to calculate potential density from potential temperature (thetao) and salinity (so)
ds["dens"] = gsw.sigma0(ds.so, ds.thetao)

ds.dens.isel(time=0).mean(dim="nlon").plot()
plt.gca().invert_yaxis()

density_grid = xcdat.create_grid(
    z=xcdat.create_axis("lev", np.linspace(6, 26, 40))
)

output = ds.regridder.vertical("so", density_grid, tool="xgcm", method="linear", target_data="dens")

output.so.isel(time=0).mean(dim="nlon").plot()
plt.gca().invert_yaxis()

2023-09-26 15:46:53,807 [WARNING]: bounds.py(add_missing_bounds:186) >> The nlat coord variable has a 'units' attribute that is not in degrees.
2023-09-26 15:46:53,807 [WARNING]: bounds.py(add_missing_bounds:186) >> The nlat coord variable has a 'units' attribute that is not in degrees.

5. Regridding using the conservative method

Related API: xarray.Dataset.regridder.vertical()

Here we will regrid the input data to the output grid using the xgcm tool and the conservative method.

We’ll transform model levels using conservative regridding. In order to perform the regridding we’ll need two grid positions, the lev coordinate is center and we”ll create the outer points using cf_xarray”s bounds_to_vertices.

ds_olev = ds.cf.bounds_to_vertices("lev").rename({"lev_vertices": "olev"})

output = ds_olev.regridder.vertical("so", output_grid, tool="xgcm", method="conservative", grid_positions={"center": "lev", "outer": "olev"})

output.so.isel(time=0).sel(lev=0, method="nearest").plot()

2023-09-26 15:47:50,128 [WARNING]: bounds.py(add_missing_bounds:186) >> The nlat coord variable has a 'units' attribute that is not in degrees.
2023-09-26 15:47:50,128 [WARNING]: bounds.py(add_missing_bounds:186) >> The nlat coord variable has a 'units' attribute that is not in degrees.

<matplotlib.collections.QuadMesh at 0x1875e1a10>

Example 2: Remapping Atmosphere Variables

1. Open dataset

# Url of data from the E3SM model in CMIP6
url_ta = 'https://esgf-data2.llnl.gov/thredds/dodsC/user_pub_work/CMIP6/CMIP/E3SM-Project/E3SM-2-0/historical/r1i1p1f1/Amon/ta/gr/v20220830/ta_Amon_E3SM-2-0_historical_r1i1p1f1_gr_185001-189912.nc'
url_cl = 'https://esgf-data2.llnl.gov/thredds/dodsC/user_pub_work/CMIP6/CMIP/E3SM-Project/E3SM-2-0/historical/r1i1p1f1/Amon/cl/gr/v20220830/cl_Amon_E3SM-2-0_historical_r1i1p1f1_gr_185001-189912.nc'

ds_ta = xcdat.open_dataset(url_ta, chunks={"time": 4}, add_bounds=["Z"])
ds_cl = xcdat.open_dataset(url_cl, chunks={"time": 4})

2. Create the output grid

Related API: xcdat.create_grid()

In this example, we will generate a grid with a linear spaced level coordinate.

output_grid = xcdat.create_grid(
    z=xcdat.create_axis("lev", np.linspace(100000, 1, 13))
)

output_grid

<xarray.Dataset>
Dimensions:   (lev: 13, bnds: 2)
Coordinates:
  * lev       (lev) float64 1e+05 9.167e+04 8.333e+04 ... 8.334e+03 1.0
Dimensions without coordinates: bnds
Data variables:
    lev_bnds  (lev, bnds) float64 1.042e+05 9.583e+04 ... 4.168e+03 -4.166e+03

xarray.Dataset

• Dimensions:
  • lev: 13
  • bnds: 2
• Coordinates: (1)
  • lev
    (lev)
    float64
    1e+05 9.167e+04 ... 8.334e+03 1.0

    axis :

    Z

    coordinate :

    vertical

    bounds :

    lev_bnds

    array([1.000000e+05, 9.166675e+04, 8.333350e+04, 7.500025e+04, 6.666700e+04,
           5.833375e+04, 5.000050e+04, 4.166725e+04, 3.333400e+04, 2.500075e+04,
           1.666750e+04, 8.334250e+03, 1.000000e+00])

• Data variables: (1)
  • lev_bnds
    (lev, bnds)
    float64
    1.042e+05 9.583e+04 ... -4.166e+03

    xcdat_bounds :

    True

    array([[104166.625,  95833.375],
           [ 95833.375,  87500.125],
           [ 87500.125,  79166.875],
           [ 79166.875,  70833.625],
           [ 70833.625,  62500.375],
           [ 62500.375,  54167.125],
           [ 54167.125,  45833.875],
           [ 45833.875,  37500.625],
           [ 37500.625,  29167.375],
           [ 29167.375,  20834.125],
           [ 20834.125,  12500.875],
           [ 12500.875,   4167.625],
           [  4167.625,  -4165.625]])

• Indexes: (1)
  • lev
    PandasIndex

    PandasIndex(Index([100000.0, 91666.75,  83333.5, 75000.25,  66667.0, 58333.75,  50000.5,
           41667.25,  33334.0, 25000.75,  16667.5,  8334.25,      1.0],
          dtype='float64', name='lev'))

• Attributes: (0)

3. Remapping air temperature on pressure levels to a set of target pressure levels.

Related API: xarray.Dataset.regridder.vertical()

Here we will regrid the input data to the output grid using the xgcm tool and the log method.

We’ll remap pressure levels.

# Remap from original pressure level to target pressure level using logarithmic interpolation
# Note: output grids can be either ascending or descending
output_ta = ds_ta.regridder.vertical("ta", output_grid, method="log")

output_ta.ta.isel(time=0, lev=0).plot()

<matplotlib.collections.QuadMesh at 0x188277b10>

4: Remap cloud fraction from model hybrid coordinate to pressure levels

Related API: xarray.Dataset.regridder.vertical()

Here we will regrid the input data to the output grid using the xgcm tool and the linear method.

We’ll remap cloud fraction into pressure space.

# Build hybrid pressure coordinate
def hybrid_coordinate(p0, a, b, ps, **kwargs):
    return a*p0 + b*ps

pressure = hybrid_coordinate(**ds_cl.data_vars)

pressure

<xarray.DataArray (lev: 72, time: 600, lat: 180, lon: 360)>
dask.array<add, shape=(72, 600, 180, 360), dtype=float64, chunksize=(72, 4, 180, 360), chunktype=numpy.ndarray>
Coordinates:
  * lev      (lev) float64 0.9985 0.9938 0.9862 ... 0.0001828 0.0001238
  * time     (time) object 1850-01-16 12:00:00 ... 1899-12-16 12:00:00
  * lat      (lat) float64 -89.5 -88.5 -87.5 -86.5 -85.5 ... 86.5 87.5 88.5 89.5
  * lon      (lon) float64 0.5 1.5 2.5 3.5 4.5 ... 355.5 356.5 357.5 358.5 359.5

xarray.DataArray

• lev: 72
• time: 600
• lat: 180
• lon: 360

• dask.array<chunksize=(72, 4, 180, 360), meta=np.ndarray>

  Array Chunk Bytes 20.86 GiB 142.38 MiB Shape (72, 600, 180, 360) (72, 4, 180, 360) Dask graph 150 chunks in 15 graph layers Data type float64 numpy.ndarray

• Coordinates: (4)
  • lev
    (lev)
    float64
    0.9985 0.9938 ... 0.0001238

    bounds :

    lev_bnds

    units :

    1

    axis :

    Z

    positive :

    down

    long_name :

    hybrid sigma pressure coordinate

    standard_name :

    atmosphere_hybrid_sigma_pressure_coordinate

    formula :

    p = a*p0 + b*ps

    formula_terms :

    p0: p0 a: a b: b ps: ps

    array([9.984964e-01, 9.937570e-01, 9.862053e-01, 9.773141e-01, 9.679111e-01,
           9.580128e-01, 9.476362e-01, 9.367990e-01, 9.255197e-01, 9.138175e-01,
           9.017118e-01, 8.892227e-01, 8.763706e-01, 8.631764e-01, 8.496612e-01,
           8.350952e-01, 8.178688e-01, 7.969527e-01, 7.723628e-01, 7.444748e-01,
           7.137046e-01, 6.804980e-01, 6.453200e-01, 6.086438e-01, 5.706249e-01,
           5.316395e-01, 4.924686e-01, 4.543958e-01, 4.188035e-01, 3.859990e-01,
           3.557641e-01, 3.278975e-01, 3.022136e-01, 2.785416e-01, 2.567237e-01,
           2.366148e-01, 2.180810e-01, 2.009989e-01, 1.852549e-01, 1.707441e-01,
           1.573699e-01, 1.450433e-01, 1.336822e-01, 1.232110e-01, 1.135600e-01,
           1.046650e-01, 9.646667e-02, 8.891054e-02, 8.194628e-02, 7.535534e-02,
           6.847639e-02, 6.101518e-02, 5.330956e-02, 4.567120e-02, 3.836628e-02,
           3.159325e-02, 2.544470e-02, 1.999634e-02, 1.533394e-02, 1.147379e-02,
           8.377404e-03, 5.968449e-03, 4.129833e-03, 2.797027e-03, 1.894352e-03,
           1.282995e-03, 8.689386e-04, 5.885091e-04, 3.985817e-04, 2.699489e-04,
           1.828292e-04, 1.238254e-04])

  • time
    (time)
    object
    1850-01-16 12:00:00 ... 1899-12-...

    bounds :

    time_bnds

    axis :

    T

    long_name :

    time

    standard_name :

    time

    _ChunkSizes :

    1

    array([cftime.DatetimeNoLeap(1850, 1, 16, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1850, 2, 15, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1850, 3, 16, 12, 0, 0, 0, has_year_zero=True),
           ...,
           cftime.DatetimeNoLeap(1899, 10, 16, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1899, 11, 16, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeNoLeap(1899, 12, 16, 12, 0, 0, 0, has_year_zero=True)],
          dtype=object)

  • lat
    (lat)
    float64
    -89.5 -88.5 -87.5 ... 88.5 89.5

    bounds :

    lat_bnds

    units :

    degrees_north

    axis :

    Y

    long_name :

    Latitude

    standard_name :

    latitude

    array([-89.5, -88.5, -87.5, -86.5, -85.5, -84.5, -83.5, -82.5, -81.5, -80.5,
           -79.5, -78.5, -77.5, -76.5, -75.5, -74.5, -73.5, -72.5, -71.5, -70.5,
           -69.5, -68.5, -67.5, -66.5, -65.5, -64.5, -63.5, -62.5, -61.5, -60.5,
           -59.5, -58.5, -57.5, -56.5, -55.5, -54.5, -53.5, -52.5, -51.5, -50.5,
           -49.5, -48.5, -47.5, -46.5, -45.5, -44.5, -43.5, -42.5, -41.5, -40.5,
           -39.5, -38.5, -37.5, -36.5, -35.5, -34.5, -33.5, -32.5, -31.5, -30.5,
           -29.5, -28.5, -27.5, -26.5, -25.5, -24.5, -23.5, -22.5, -21.5, -20.5,
           -19.5, -18.5, -17.5, -16.5, -15.5, -14.5, -13.5, -12.5, -11.5, -10.5,
            -9.5,  -8.5,  -7.5,  -6.5,  -5.5,  -4.5,  -3.5,  -2.5,  -1.5,  -0.5,
             0.5,   1.5,   2.5,   3.5,   4.5,   5.5,   6.5,   7.5,   8.5,   9.5,
            10.5,  11.5,  12.5,  13.5,  14.5,  15.5,  16.5,  17.5,  18.5,  19.5,
            20.5,  21.5,  22.5,  23.5,  24.5,  25.5,  26.5,  27.5,  28.5,  29.5,
            30.5,  31.5,  32.5,  33.5,  34.5,  35.5,  36.5,  37.5,  38.5,  39.5,
            40.5,  41.5,  42.5,  43.5,  44.5,  45.5,  46.5,  47.5,  48.5,  49.5,
            50.5,  51.5,  52.5,  53.5,  54.5,  55.5,  56.5,  57.5,  58.5,  59.5,
            60.5,  61.5,  62.5,  63.5,  64.5,  65.5,  66.5,  67.5,  68.5,  69.5,
            70.5,  71.5,  72.5,  73.5,  74.5,  75.5,  76.5,  77.5,  78.5,  79.5,
            80.5,  81.5,  82.5,  83.5,  84.5,  85.5,  86.5,  87.5,  88.5,  89.5])

  • lon
    (lon)
    float64
    0.5 1.5 2.5 ... 357.5 358.5 359.5

    bounds :

    lon_bnds

    units :

    degrees_east

    axis :

    X

    long_name :

    Longitude

    standard_name :

    longitude

    array([  0.5,   1.5,   2.5, ..., 357.5, 358.5, 359.5])

• Indexes: (4)
  • lev
    PandasIndex

    PandasIndex(Index([    0.9984964394917621,     0.9937569563042221,     0.9862053331610676,
                0.977314063869485,     0.9679110804007753,     0.9580127663713389,
               0.9476361588543765,     0.9367989948227038,      0.925519748147335,
                0.913817504420407,     0.9017118332902198,     0.8892227360962712,
               0.8763706483341046,       0.86317642564934,     0.8496612491105949,
               0.8350951708866112,     0.8178688234772361,      0.796952749718736,
               0.7723628379433543,     0.7444748293231951,     0.7137046383204335,
               0.6804980444382399,     0.6453199679615715,     0.6086437744215841,
               0.5706249031033329,      0.531639531357463,    0.49246857402252703,
               0.4543958139891781,     0.4188034628792632,     0.3859990239199109,
              0.35576412943838015,    0.32789750608917023,     0.3022136401218828,
               0.2785415583098347,     0.2567236859052588,    0.23661478658749036,
              0.21808099703300082,    0.20099893987110398,    0.18525490222915122,
              0.17074407952158202,    0.15736987784007928,    0.14504326553270042,
               0.1336821799321061,      0.123210991427466,    0.11356000142146327,
              0.10466496728675462,    0.09646667344493705,    0.08891054314956091,
              0.08194627512473471,    0.07535533589726098,    0.06847639022971255,
             0.061015181670369456,   0.053309561434392534,   0.045671197939426804,
              0.03836628309442701,    0.03159325129128525,   0.025444696509961597,
             0.019996337978167075,   0.015333938222216858,   0.011473787230534558,
              0.00837740437894177,   0.005968449367118077,   0.004129833155024836,
            0.0027970269352932593,  0.0018943524794063876,   0.001282994908254925,
            0.0008689385700405001,   0.000588509146564753, 0.00039858170362288225,
           0.00026994886212093746, 0.00018282923550684864, 0.00012382541305561677],
          dtype='float64', name='lev'))

  • time
    PandasIndex

    PandasIndex(CFTimeIndex([1850-01-16 12:00:00, 1850-02-15 00:00:00, 1850-03-16 12:00:00,
                 1850-04-16 00:00:00, 1850-05-16 12:00:00, 1850-06-16 00:00:00,
                 1850-07-16 12:00:00, 1850-08-16 12:00:00, 1850-09-16 00:00:00,
                 1850-10-16 12:00:00,
                 ...
                 1899-03-16 12:00:00, 1899-04-16 00:00:00, 1899-05-16 12:00:00,
                 1899-06-16 00:00:00, 1899-07-16 12:00:00, 1899-08-16 12:00:00,
                 1899-09-16 00:00:00, 1899-10-16 12:00:00, 1899-11-16 00:00:00,
                 1899-12-16 12:00:00],
                dtype='object', length=600, calendar='noleap', freq='None'))

  • lat
    PandasIndex

    PandasIndex(Index([-89.5, -88.5, -87.5, -86.5, -85.5, -84.5, -83.5, -82.5, -81.5, -80.5,
           ...
            80.5,  81.5,  82.5,  83.5,  84.5,  85.5,  86.5,  87.5,  88.5,  89.5],
          dtype='float64', name='lat', length=180))

  • lon
    PandasIndex

    PandasIndex(Index([  0.5,   1.5,   2.5,   3.5,   4.5,   5.5,   6.5,   7.5,   8.5,   9.5,
           ...
           350.5, 351.5, 352.5, 353.5, 354.5, 355.5, 356.5, 357.5, 358.5, 359.5],
          dtype='float64', name='lon', length=360))

• Attributes: (0)

new_pressure_grid = xcdat.create_grid(
    z=xcdat.create_axis("lev", np.linspace(100000, 1, 13))
)

output_cl = ds_cl.regridder.vertical("cl", new_pressure_grid, method="linear", target_data=pressure)

output_cl.cl.isel(time=0, lev=0).plot()

<matplotlib.collections.QuadMesh at 0x19974eed0>

output_cl.cl.isel(time=0).mean(dim='lon').plot()
plt.gca().invert_yaxis()