history files usage, case metadata generated from CESM CASEROOT#

import ast
import os
import os.path
import pprint
import warnings

from dask_jobqueue import PBSCluster
from distributed import Client
from distributed.worker import get_client
import intake

from esm_catalog_utils import query_from_caseroot, caseroot_to_esm_datastore
%load_ext autoreload
%autoreload 2

Specification of Input#

caseroot = "/glade/work/klindsay/cesm20_cases/B1850/b.e21.B1850.f19_g17.ex_output.001"

Obtain Computational Resources#

user = os.getenv('USER')
logdir = f'/glade/scratch/{user}/dask_tmp'
with warnings.catch_warnings():
    warnings.filterwarnings(action="ignore", message=".*already in use", module=".*node")
    cluster = PBSCluster(
        cores = 1,
        memory = '4GiB',
        processes = 1,
        local_directory = logdir,
        log_directory = logdir,
        resource_spec = 'select=1:ncpus=1:mem=4GB',
        queue = 'casper',
        account = 'P93300070',
        walltime = '0:30:00',
        interface = 'ib0',
    )

cluster.scale(12)

client = Client(cluster)
client

Client

Client-d04fd975-f4f1-11ed-95d1-3cecef1aca66

Connection method: Cluster object Cluster type: dask_jobqueue.PBSCluster
Dashboard: https://jupyterhub.hpc.ucar.edu/stable/user/klindsay/esm_catalog_utils/proxy/8787/status

Cluster Info

Generate and save catalog#

%%time
print("creating esm_datastore")
esm_datastore = caseroot_to_esm_datastore(caseroot, use_dask=True)
print(esm_datastore)
creating esm_datastore
generating esm_datastore for b.e21.B1850.f19_g17.ex_output.001
<sample catalog with 13 dataset(s) from 351 asset(s)>
CPU times: user 1.18 s, sys: 248 ms, total: 1.42 s
Wall time: 12.7 s
df = esm_datastore.df
vals = df.frequency.unique() 
print(vals)
['month_1' 'day_1' 'year_1' 'unknown' 'once']
case = query_from_caseroot(caseroot, "CASE")
datastore_directory = os.path.join(os.getcwd(), "generated")
os.makedirs(datastore_directory, exist_ok=True)
esm_datastore.serialize(name=case, directory=datastore_directory, catalog_type="file")
Successfully wrote ESM catalog json file to: file:///glade/work/klindsay/analysis/esm_catalog_utils/notebooks/generated/b.e21.B1850.f19_g17.ex_output.001.json

Generate catalog, in update mode#

%%time
path = os.path.join(datastore_directory, f"{case}.json")
print("updating esm_datastore")
esm_datastore = intake.open_esm_datastore(path)
print(esm_datastore)
esm_datastore = caseroot_to_esm_datastore(
    caseroot, esm_datastore_in=esm_datastore, use_dask=True
)
print(esm_datastore)
updating esm_datastore
<b.e21.B1850.f19_g17.ex_output.001 catalog with 13 dataset(s) from 351 asset(s)>
appending esm_datastore for b.e21.B1850.f19_g17.ex_output.001
<b.e21.B1850.f19_g17.ex_output.001 catalog with 13 dataset(s) from 351 asset(s)>
CPU times: user 537 ms, sys: 29.2 ms, total: 567 ms
Wall time: 4.53 s

Load catalog#

use read_csv_kwargs argument to enable parsing of list output in varname column#

esm_datastore = intake.open_esm_datastore(
    path, read_csv_kwargs={"converters": {"varname": ast.literal_eval}}
)
print(esm_datastore)
esm_datastore.df.columns
<b.e21.B1850.f19_g17.ex_output.001 catalog with 13 dataset(s) from 351 asset(s)>
Index(['case', 'scomp', 'component', 'path', 'stream', 'datestring',
       'frequency', 'date_start', 'date_end', 'varname', 'size'],
      dtype='object')
df = esm_datastore.df
vals = df.frequency.unique() 
print(vals)
['month_1' 'day_1' 'year_1' 'unknown' 'once']

verify that dataset_dict can be created for each scomp#

## "coords": "minimal" required for CLM history files if output spans multiple submissions
## "compat": "override" improves performance of combine step
##   and requires "coords": "minimal"
xarray_combine_by_coords_kwargs = {
    "coords": "minimal", "compat": "override", "data_vars": "minimal"
}
for scomp in df.scomp.unique():
    print(f"scomp={scomp}")
    subset = esm_datastore.search(scomp=scomp)
    print(subset)

    with warnings.catch_warnings():
        warnings.filterwarnings(
            action="ignore", message=".*single-machine", module=".*base"
        )
        # ds_dict = subset.to_dataset_dict(xarray_open_kwargs={"drop_variables": ["average_T1", "average_T2"]})
        ds_dict = subset.to_dataset_dict(
            xarray_combine_by_coords_kwargs=xarray_combine_by_coords_kwargs
        )
        pprint.pprint(list(ds_dict.keys()))
        key = list(ds_dict)[0]
        ds = ds_dict[key]
        print(ds)

    print("============================================================")
scomp=cam
<b.e21.B1850.f19_g17.ex_output.001 catalog with 2 dataset(s) from 52 asset(s)>

--> The keys in the returned dictionary of datasets are constructed as follows:
	'case.scomp.component.stream.frequency'
100.00% [2/2 00:12<00:00]
['b.e21.B1850.f19_g17.ex_output.001.cam.atm.i.day_1',
 'b.e21.B1850.f19_g17.ex_output.001.cam.atm.h0.month_1']
<xarray.Dataset>
Dimensions:       (slat: 95, lon: 144, lat: 96, slon: 144, lev: 32, ilev: 33,
                   time: 4, nbnd: 2)
Coordinates: (12/20)
  * slat          (slat) float64 -89.05 -87.16 -85.26 ... 85.26 87.16 89.05
  * lon           (lon) float64 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
    w_stag        (slat) float64 dask.array<chunksize=(95,), meta=np.ndarray>
  * lat           (lat) float64 -90.0 -88.11 -86.21 -84.32 ... 86.21 88.11 90.0
  * slon          (slon) float64 -1.25 1.25 3.75 6.25 ... 351.2 353.8 356.2
    gw            (lat) float64 dask.array<chunksize=(96,), meta=np.ndarray>
    ...            ...
    time_bnds     (time, nbnd) object dask.array<chunksize=(1, 2), meta=np.ndarray>
    ndbase        int32 ...
    nsbase        int32 ...
    nbdate        int32 ...
    nbsec         int32 ...
    mdt           int32 ...
Dimensions without coordinates: nbnd
Data variables: (12/48)
    date          (time) int32 dask.array<chunksize=(1,), meta=np.ndarray>
    datesec       (time) int32 dask.array<chunksize=(1,), meta=np.ndarray>
    date_written  (time) |S8 dask.array<chunksize=(1,), meta=np.ndarray>
    time_written  (time) |S8 dask.array<chunksize=(1,), meta=np.ndarray>
    ndcur         (time) int32 dask.array<chunksize=(1,), meta=np.ndarray>
    nscur         (time) int32 dask.array<chunksize=(1,), meta=np.ndarray>
    ...            ...
    pom_a4        (time, lev, lat, lon) float64 dask.array<chunksize=(1, 32, 96, 144), meta=np.ndarray>
    so4_a1        (time, lev, lat, lon) float64 dask.array<chunksize=(1, 32, 96, 144), meta=np.ndarray>
    so4_a2        (time, lev, lat, lon) float64 dask.array<chunksize=(1, 32, 96, 144), meta=np.ndarray>
    so4_a3        (time, lev, lat, lon) float64 dask.array<chunksize=(1, 32, 96, 144), meta=np.ndarray>
    soa_a1        (time, lev, lat, lon) float64 dask.array<chunksize=(1, 32, 96, 144), meta=np.ndarray>
    soa_a2        (time, lev, lat, lon) float64 dask.array<chunksize=(1, 32, 96, 144), meta=np.ndarray>
Attributes: (12/17)
    Conventions:                     CF-1.0
    source:                          CAM
    case:                            b.e21.B1850.f19_g17.ex_output.001
    logname:                         klindsay
    initial_file:                    b.e20.B1850.f19_g17.release_cesm2_1_0.02...
    topography_file:                 /glade/p/cesmdata/cseg/inputdata/atm/cam...
    ...                              ...
    intake_esm_attrs:component:      atm
    intake_esm_attrs:stream:         i
    intake_esm_attrs:frequency:      day_1
    intake_esm_attrs:varname:        date,datesec,date_written,time_written,n...
    intake_esm_attrs:_data_format_:  netcdf
    intake_esm_dataset_key:          b.e21.B1850.f19_g17.ex_output.001.cam.at...
============================================================
scomp=cism
<b.e21.B1850.f19_g17.ex_output.001 catalog with 2 dataset(s) from 5 asset(s)>

--> The keys in the returned dictionary of datasets are constructed as follows:
	'case.scomp.component.stream.frequency'
100.00% [2/2 00:03<00:00]
['b.e21.B1850.f19_g17.ex_output.001.cism.glc.initial_hist.unknown',
 'b.e21.B1850.f19_g17.ex_output.001.cism.glc.h.year_1']
<xarray.Dataset>
Dimensions:        (time: 1, level: 11, lithoz: 20, staglevel: 10,
                    stagwbndlevel: 12, x0: 415, x1: 416, y0: 703, y1: 704)
Coordinates:
  * time           (time) object 0001-01-01 00:00:00
  * level          (level) float64 0.0 0.2314 0.4074 0.5444 ... 0.9218 0.964 1.0
  * lithoz         (lithoz) float64 9.969e+36 9.969e+36 ... 9.969e+36 9.969e+36
  * staglevel      (staglevel) float64 0.1157 0.3194 0.4759 ... 0.9429 0.982
  * stagwbndlevel  (stagwbndlevel) float64 0.0 0.1157 0.3194 ... 0.982 1.0
  * x0             (x0) float64 -7.11e+05 -7.07e+05 ... 9.41e+05 9.45e+05
  * x1             (x1) float64 -7.13e+05 -7.09e+05 ... 9.43e+05 9.47e+05
  * y0             (y0) float64 -3.394e+06 -3.39e+06 ... -5.9e+05 -5.86e+05
  * y1             (y1) float64 -3.396e+06 -3.392e+06 ... -5.88e+05 -5.84e+05
Data variables:
    internal_time  (time) object dask.array<chunksize=(1,), meta=np.ndarray>
    tstep_count    (time) int32 dask.array<chunksize=(1,), meta=np.ndarray>
    mapping        |S1 ...
    artm           (time, y1, x1) float64 dask.array<chunksize=(1, 704, 416), meta=np.ndarray>
    smb            (time, y1, x1) float64 dask.array<chunksize=(1, 704, 416), meta=np.ndarray>
    thk            (time, y1, x1) float64 dask.array<chunksize=(1, 704, 416), meta=np.ndarray>
    topg           (time, y1, x1) float64 dask.array<chunksize=(1, 704, 416), meta=np.ndarray>
    usurf          (time, y1, x1) float64 dask.array<chunksize=(1, 704, 416), meta=np.ndarray>
Attributes: (12/28)
    Conventions:                     CF-1.3
    title:                           
    institution:                     
    source:                          
    history:                         
    references:                      
    ...                              ...
    intake_esm_attrs:date_start:     0001-01-01
    intake_esm_attrs:date_end:       0001-01-01
    intake_esm_attrs:varname:        internal_time,tstep_count,artm,smb,thk,t...
    intake_esm_attrs:size:           11736996
    intake_esm_attrs:_data_format_:  netcdf
    intake_esm_dataset_key:          b.e21.B1850.f19_g17.ex_output.001.cism.g...
============================================================
scomp=cice
<b.e21.B1850.f19_g17.ex_output.001 catalog with 1 dataset(s) from 48 asset(s)>

--> The keys in the returned dictionary of datasets are constructed as follows:
	'case.scomp.component.stream.frequency'
100.00% [1/1 00:02<00:00]
['b.e21.B1850.f19_g17.ex_output.001.cice.ice.h.month_1']
<xarray.Dataset>
Dimensions:      (time: 48, d2: 2, nj: 384, ni: 320, nc: 5, nkice: 8,
                  nksnow: 3, nkbio: 5, nvertices: 4)
Coordinates: (12/25)
  * time         (time) object 0001-02-01 00:00:00 ... 0005-01-01 00:00:00
    time_bounds  (time, d2) object dask.array<chunksize=(1, 2), meta=np.ndarray>
    TLON         (nj, ni) float32 dask.array<chunksize=(384, 320), meta=np.ndarray>
    TLAT         (nj, ni) float32 dask.array<chunksize=(384, 320), meta=np.ndarray>
    ULON         (nj, ni) float32 dask.array<chunksize=(384, 320), meta=np.ndarray>
    ULAT         (nj, ni) float32 dask.array<chunksize=(384, 320), meta=np.ndarray>
    ...           ...
    ANGLE        (nj, ni) float32 dask.array<chunksize=(384, 320), meta=np.ndarray>
    ANGLET       (nj, ni) float32 dask.array<chunksize=(384, 320), meta=np.ndarray>
    lont_bounds  (nj, ni, nvertices) float32 dask.array<chunksize=(384, 320, 4), meta=np.ndarray>
    latt_bounds  (nj, ni, nvertices) float32 dask.array<chunksize=(384, 320, 4), meta=np.ndarray>
    lonu_bounds  (nj, ni, nvertices) float32 dask.array<chunksize=(384, 320, 4), meta=np.ndarray>
    latu_bounds  (nj, ni, nvertices) float32 dask.array<chunksize=(384, 320, 4), meta=np.ndarray>
Dimensions without coordinates: d2, nj, ni, nc, nkice, nksnow, nkbio, nvertices
Data variables: (12/68)
    hi           (time, nj, ni) float32 dask.array<chunksize=(1, 384, 320), meta=np.ndarray>
    hs           (time, nj, ni) float32 dask.array<chunksize=(1, 384, 320), meta=np.ndarray>
    snowfrac     (time, nj, ni) float32 dask.array<chunksize=(1, 384, 320), meta=np.ndarray>
    Tsfc         (time, nj, ni) float32 dask.array<chunksize=(1, 384, 320), meta=np.ndarray>
    aice         (time, nj, ni) float32 dask.array<chunksize=(1, 384, 320), meta=np.ndarray>
    uvel         (time, nj, ni) float32 dask.array<chunksize=(1, 384, 320), meta=np.ndarray>
    ...           ...
    aicen        (time, nc, nj, ni) float32 dask.array<chunksize=(1, 5, 384, 320), meta=np.ndarray>
    vicen        (time, nc, nj, ni) float32 dask.array<chunksize=(1, 5, 384, 320), meta=np.ndarray>
    vsnon        (time, nc, nj, ni) float32 dask.array<chunksize=(1, 5, 384, 320), meta=np.ndarray>
    fswsfcn      (time, nc, nj, ni) float32 dask.array<chunksize=(1, 5, 384, 320), meta=np.ndarray>
    fswintn      (time, nc, nj, ni) float32 dask.array<chunksize=(1, 5, 384, 320), meta=np.ndarray>
    fswthrun     (time, nc, nj, ni) float32 dask.array<chunksize=(1, 5, 384, 320), meta=np.ndarray>
Attributes: (12/19)
    title:                           b.e21.B1850.f19_g17.ex_output.001
    contents:                        Diagnostic and Prognostic Variables
    source:                          Los Alamos Sea Ice Model (CICE) Version 5
    time_period_freq:                month_1
    model_doi_url:                   https://doi.org/10.5065/D67H1H0V
    comment:                         All years have exactly 365 days
    ...                              ...
    intake_esm_attrs:stream:         h
    intake_esm_attrs:frequency:      month_1
    intake_esm_attrs:varname:        hi,hs,snowfrac,Tsfc,aice,uvel,vvel,uatm,...
    intake_esm_attrs:size:           60485728
    intake_esm_attrs:_data_format_:  netcdf
    intake_esm_dataset_key:          b.e21.B1850.f19_g17.ex_output.001.cice.i...
============================================================
scomp=clm2
<b.e21.B1850.f19_g17.ex_output.001 catalog with 1 dataset(s) from 48 asset(s)>

--> The keys in the returned dictionary of datasets are constructed as follows:
	'case.scomp.component.stream.frequency'
100.00% [1/1 00:14<00:00]
['b.e21.B1850.f19_g17.ex_output.001.clm2.lnd.h0.month_1']
<xarray.Dataset>
Dimensions:                          (levgrnd: 25, levlak: 10, levdcmp: 25,
                                      time: 48, hist_interval: 2, lon: 144,
                                      lat: 96, levsoi: 20, cft: 64,
                                      glc_nec: 10, ltype: 9, natpft: 15,
                                      nvegwcs: 4)
Coordinates: (12/20)
  * levgrnd                          (levgrnd) float32 0.01 0.04 ... 28.87 42.0
  * levlak                           (levlak) float32 0.05 0.6 ... 34.33 44.78
  * levdcmp                          (levdcmp) float32 0.01 0.04 ... 28.87 42.0
  * time                             (time) object 0001-02-01 00:00:00 ... 00...
    time_bounds                      (time, hist_interval) object dask.array<chunksize=(1, 2), meta=np.ndarray>
  * lon                              (lon) float32 0.0 2.5 5.0 ... 355.0 357.5
    ...                               ...
    WATSAT                           (levgrnd, lat, lon) float32 dask.array<chunksize=(25, 96, 144), meta=np.ndarray>
    SUCSAT                           (levgrnd, lat, lon) float32 dask.array<chunksize=(25, 96, 144), meta=np.ndarray>
    BSW                              (levgrnd, lat, lon) float32 dask.array<chunksize=(25, 96, 144), meta=np.ndarray>
    HKSAT                            (levgrnd, lat, lon) float32 dask.array<chunksize=(25, 96, 144), meta=np.ndarray>
    ZLAKE                            (levlak, lat, lon) float32 dask.array<chunksize=(10, 96, 144), meta=np.ndarray>
    DZLAKE                           (levlak, lat, lon) float32 dask.array<chunksize=(10, 96, 144), meta=np.ndarray>
Dimensions without coordinates: hist_interval, levsoi, cft, glc_nec, ltype,
                                natpft, nvegwcs
Data variables: (12/465)
    mcdate                           (time) int32 dask.array<chunksize=(1,), meta=np.ndarray>
    mcsec                            (time) int32 dask.array<chunksize=(1,), meta=np.ndarray>
    mdcur                            (time) int32 dask.array<chunksize=(1,), meta=np.ndarray>
    mscur                            (time) int32 dask.array<chunksize=(1,), meta=np.ndarray>
    nstep                            (time) int32 dask.array<chunksize=(1,), meta=np.ndarray>
    date_written                     (time) |S16 dask.array<chunksize=(1,), meta=np.ndarray>
    ...                               ...
    XSMRPOOL                         (time, lat, lon) float32 dask.array<chunksize=(1, 96, 144), meta=np.ndarray>
    XSMRPOOL_RECOVER                 (time, lat, lon) float32 dask.array<chunksize=(1, 96, 144), meta=np.ndarray>
    ZBOT                             (time, lat, lon) float32 dask.array<chunksize=(1, 96, 144), meta=np.ndarray>
    ZWT                              (time, lat, lon) float32 dask.array<chunksize=(1, 96, 144), meta=np.ndarray>
    ZWT_CH4_UNSAT                    (time, lat, lon) float32 dask.array<chunksize=(1, 96, 144), meta=np.ndarray>
    ZWT_PERCH                        (time, lat, lon) float32 dask.array<chunksize=(1, 96, 144), meta=np.ndarray>
Attributes: (12/108)
    title:                                     CLM History file information
    comment:                                   NOTE: None of the variables ar...
    Conventions:                               CF-1.0
    source:                                    Community Land Model CLM4.0
    hostname:                                  cheyenne
    username:                                  klindsay
    ...                                        ...
    intake_esm_attrs:stream:                   h0
    intake_esm_attrs:frequency:                month_1
    intake_esm_attrs:varname:                  mcdate,mcsec,mdcur,mscur,nstep...
    intake_esm_attrs:_data_format_:            netcdf
    Time_constant_3Dvars:                      ZSOI:DZSOI:WATSAT:SUCSAT:BSW:H...
    intake_esm_dataset_key:                    b.e21.B1850.f19_g17.ex_output....
============================================================
scomp=pop
<b.e21.B1850.f19_g17.ex_output.001 catalog with 6 dataset(s) from 150 asset(s)>

--> The keys in the returned dictionary of datasets are constructed as follows:
	'case.scomp.component.stream.frequency'
50.00% [3/6 00:04<00:04]

plot timeseries of daily FG_CO2 from the ocean model#

varname = 'FG_CO2_2'
subset = esm_datastore.search(frequency='day_1', varname=varname)
print(subset)
%%time
with warnings.catch_warnings():
    warnings.filterwarnings(
        action="ignore", message=".*single-machine", module=".*base"
    )
    ds_dict = subset.to_dataset_dict(
        xarray_combine_by_coords_kwargs=xarray_combine_by_coords_kwargs
    )
key = list(ds_dict)[0]
ds = ds_dict[key]
ds
dims = ds[varname].dims[-2:]
ds[varname].isel(time=slice(0,365*3)).mean(dims).plot()

Release Computational Resources#

client.close()
cluster.close()