Comparing Atmospheric Model Output with Observations Using Intake-ESM#

Comparing models and observations is a critical component of climate diagnostic packages. This process can be challenging though - given the number of observational datasets to compare against, and the difference in spatiotemporal resolutions. In the previous iteration of the diagnostics package used for atmospheric data from the Community Earth System Model (CESM), they used pre-computed, observational datasets stored in a directory on the GLADE filesystem (/glade/p/cesm/amwg/amwg_diagnostics/obs_data)

Within this example, we walk though generating an intake-esm catalog from the observational data, reading in CESM data, and compare models and observations

Imports#

import ast
import glob
import pathlib
import traceback
import warnings
from datetime import datetime

import cartopy.crs as ccrs
import geocat.comp
import geoviews.feature as gf
import holoviews as hv
import hvplot
import hvplot.xarray
import intake
import matplotlib.pyplot as plt
import pandas as pd
import xarray as xr
from distributed import Client
from ecgtools import Builder
from ecgtools.builder import INVALID_ASSET, TRACEBACK
from ncar_jobqueue import NCARCluster

hv.extension('bokeh')

Build an `Intake-ESM` catalog from the observational dataset#

Before we jump into building the catalog, let’s take a look at what data we are working with

Build a parser using `ecgtools`#

We can use the ecgtools package to help build an intake-esm catalog for this dataset!

We start first by using pathlib to generically look at the filepaths

path = pathlib.Path(files[0])
path.stem

'AIRS_01_climo'

This path can be split using .split('_'), separates the path into the following:

observational dataset source
month number or season
“climo”

path.stem.split('_')

['AIRS', '01', 'climo']

We can also gather useful insight by opening the file!

ds = xr.open_dataset(files[0])

Let’s look at the variable “Temperature” (T)

You’ll notice that we have additional attributes about the dataset which we can use, including:

units
long_name
climatology, which includes more information about what dates were used

Adding the parsing function#

def parse_amwg_obs(file):
    """Atmospheric observational data used within the AMWG diagnostic package

    Parameters
    ----------
    file: str
      filepath from the data directory including observational climatologies

    Returns
    -------
    info: dict
       Dictionary with information to add as columns in the data catalog

    """
    file = pathlib.Path(file)
    info_list = []
    info = {}

    try:
        stem = file.stem
        split = stem.split('_')
        info['source'] = split[0]
        temporal = split[-2]
        if len(split[-2]) == 2:
            month_number = int(temporal)
            info['time_period'] = 'monthly'
            info['temporal'] = datetime(2020, month_number, 1).strftime('%b').upper()

        elif len(temporal) == 3:
            info['time_period'] = 'seasonal'
            info['temporal'] = temporal

        else:
            info['time_period'] = np.nan
            info['temporal'] = np.nan

        with xr.open_dataset(file, chunks={}, decode_times=False) as ds:
            variables = [v for v, da in ds.variables.items() if 'time' in da.dims]
            info['variables'] = variables

            info['path'] = str(path)

        return info

    except Exception:
        return {INVALID_ASSET: file, TRACEBACK: traceback.format_exc()}

'time' in ds['T'].dims

True

info = parse_amwg_obs(files[0])
info

{'source': 'AIRS',
 'time_period': 'monthly',
 'temporal': 'JAN',
 'variables': ['T', 'time', 'RELHUM', 'O3', 'SHUM'],
 'path': '/glade/p/cesm/amwg/amwg_diagnostics/obs_data/AIRS_01_climo.nc'}

pd.DataFrame(info)

	source	time_period	temporal	variables	path
0	AIRS	monthly	JAN	T	/glade/p/cesm/amwg/amwg_diagnostics/obs_data/A...
1	AIRS	monthly	JAN	time	/glade/p/cesm/amwg/amwg_diagnostics/obs_data/A...
2	AIRS	monthly	JAN	RELHUM	/glade/p/cesm/amwg/amwg_diagnostics/obs_data/A...
3	AIRS	monthly	JAN	O3	/glade/p/cesm/amwg/amwg_diagnostics/obs_data/A...
4	AIRS	monthly	JAN	SHUM	/glade/p/cesm/amwg/amwg_diagnostics/obs_data/A...

We can test this on a single file!

This dictionary can easily be converted into a pandas.Dataframe

pd.DataFrame([info, info])

	source	time_period	temporal	variables	path
0	AIRS	monthly	JAN	[T, time, RELHUM, O3, SHUM]	/glade/p/cesm/amwg/amwg_diagnostics/obs_data/A...
1	AIRS	monthly	JAN	[T, time, RELHUM, O3, SHUM]	/glade/p/cesm/amwg/amwg_diagnostics/obs_data/A...

Use `ecgtools` to build the catalog#

Now that we have our parser ready to go, we can use ecgtools to build the catalog

b = Builder(
    # Directory with the output
    f"/glade/p/cesm/amwg/amwg_diagnostics/obs_data/",
    depth=3,
    # Number of jobs to execute - should be equal to # threads you are using
    njobs=-1,
)

Once we set up the builder, can call .build, passing in the parser we built

b.build(parse_amwg_obs)

b.df

	source	time_period	temporal	variables	path
0	AIRS	monthly	JAN	[T, time, RELHUM, O3, SHUM]	/glade/p/cesm/amwg/amwg_diagnostics/obs_data/A...
1	AIRS	monthly	FEB	[T, time, RELHUM, O3, SHUM]	/glade/p/cesm/amwg/amwg_diagnostics/obs_data/A...
2	AIRS	monthly	MAR	[T, time, RELHUM, O3, SHUM]	/glade/p/cesm/amwg/amwg_diagnostics/obs_data/A...
3	AIRS	monthly	APR	[T, time, RELHUM, O3, SHUM]	/glade/p/cesm/amwg/amwg_diagnostics/obs_data/A...
4	AIRS	monthly	MAY	[T, time, RELHUM, O3, SHUM]	/glade/p/cesm/amwg/amwg_diagnostics/obs_data/A...
...	...	...	...	...	...
3778	mlsw	seasonal	ANN	[T, time, Z3, H2O, O3, RELHUM]	/glade/p/cesm/amwg/amwg_diagnostics/obs_data/A...
3779	mlsw	seasonal	DJF	[T, time, Z3, H2O, O3, RELHUM]	/glade/p/cesm/amwg/amwg_diagnostics/obs_data/A...
3780	mlsw	seasonal	JJA	[T, time, Z3, H2O, O3, RELHUM]	/glade/p/cesm/amwg/amwg_diagnostics/obs_data/A...
3781	mlsw	seasonal	MAM	[T, time, Z3, H2O, O3, RELHUM]	/glade/p/cesm/amwg/amwg_diagnostics/obs_data/A...
3782	mlsw	seasonal	SON	[T, time, Z3, H2O, O3, RELHUM]	/glade/p/cesm/amwg/amwg_diagnostics/obs_data/A...

1243 rows × 5 columns

Save the catalog#

b.save(
    # File path - could save as .csv (uncompressed csv) or .csv.gz (compressed csv)
    "/glade/work/mgrover/intake-esm-catalogs/amwg_obs_datasets.csv",
    # Column name including filepath
    path_column_name='path',
    # Column name including variables
    variable_column_name='variables',
    # Data file format - could be netcdf or zarr (in this case, netcdf)
    data_format="netcdf",
    # Which attributes to groupby when reading in variables using intake-esm
    groupby_attrs=["source", "time_period"],
    # Aggregations which are fed into xarray when reading in data using intake
    aggregations=[
        {
            'type': 'join_new',
            'attribute_name': 'temporal',
            'options': {'coords': 'minimal', 'compat': 'override'},
        },
    ],
)

Saved catalog location: /glade/work/mgrover/intake-esm-catalogs/amwg_obs_datasets.json and /glade/work/mgrover/intake-esm-catalogs/amwg_obs_datasets.csv

/glade/scratch/mgrover/ipykernel_140729/3449059562.py:1: UserWarning: Unable to parse 2541 assets/files. A list of these assets can be found in /glade/work/mgrover/intake-esm-catalogs/invalid_assets_amwg_obs_datasets.csv.
  b.save(

Read in Observational Data from the Catalog#

We use intake-esm to read in the observational data from the catalog we just created

obs_catalog = intake.open_esm_datastore(
    "/glade/work/mgrover/intake-esm-catalogs/amwg_obs_datasets.json",
    csv_kwargs={"converters": {"variables": ast.literal_eval}},
    sep="/",
)

We are interested in Temperature (T), from the AIRS dataset, which comes from NASA

obs_catalog_subset = obs_catalog.search(variables='T', source='AIRS')

dsets = obs_catalog_subset.to_dataset_dict()

--> The keys in the returned dictionary of datasets are constructed as follows:
	'source/time_period'

100.00% [2/2 00:00<00:00]

Our dictionary of datasets has two keys, one for seasonal and the other for monthly climatologies

dsets.keys()

dict_keys(['AIRS/seasonal', 'AIRS/monthly'])

seasonal_obs_ds = dsets['AIRS/seasonal']
monthly_obs_ds = dsets['AIRS/monthly']

Plot Observational Temperature Climatologies#

Now that we read in our data, we can plot up our results!

Plot Seasonal Temperature Climatologies from Observations#

seasonal_obs_plot = (
    seasonal_obs_ds.isel(time=0, lev=range(3)).T.hvplot.quadmesh(
        rasterize=True, groupby=['lev', 'temporal'], cmap='magma', projection=ccrs.Robinson()
    )
    * gf.coastline
)
seasonal_obs_plot

Plot Monthly Temperature Climatologies from Observations#

monthly_obs_plot = (
    monthly_obs_ds.isel(time=0, lev=range(3)).T.hvplot.quadmesh(
        rasterize=True, groupby=['lev', 'temporal'], cmap='magma', projection=ccrs.Robinson()
    )
    * gf.coastline
)
monthly_obs_plot

27 August 2021

Recent Posts

Archives

Comparing Atmospheric Model Output with Observations Using Intake-ESM#

Imports#

Build an `Intake-ESM` catalog from the observational dataset#

Investigate the File Structure and Sample Dataset#

Build a parser using `ecgtools`#

Adding the parsing function#

Use `ecgtools` to build the catalog#

Save the catalog#

Read in Observational Data from the Catalog#

Plot Observational Temperature Climatologies#

Plot Seasonal Temperature Climatologies from Observations#

Plot Monthly Temperature Climatologies from Observations#

27 August 2021

Recent Posts

Archives

Comparing Atmospheric Model Output with Observations Using Intake-ESM#

Imports#

Build an Intake-ESM catalog from the observational dataset#

Investigate the File Structure and Sample Dataset#

Build a parser using ecgtools#

Adding the parsing function#

Use ecgtools to build the catalog#

Save the catalog#

Read in Observational Data from the Catalog#

Plot Observational Temperature Climatologies#

Plot Seasonal Temperature Climatologies from Observations#

Plot Monthly Temperature Climatologies from Observations#

Build an `Intake-ESM` catalog from the observational dataset#

Build a parser using `ecgtools`#

Use `ecgtools` to build the catalog#