Advanced example: define advanced metrics

%%javascript
// leave this in to disable autoscroll in Jupyter notebook
IPython.OutputArea.prototype._should_scroll = function(lines) {
    return false;
}

Imports

import fcpy

Load data

# the `load_dataset' option is currently unavalable due to data handling system move.
# ds = fcpy.load_dataset(model='atmospheric-model', var_names=['2t', 'q', 'sst'], levels=list(range(1, 10)))

# The `data` folder contains two sample NetCDF files.
# Here we load specific humidity from CAMS at 32 bits
fpath = "../data/cams_q_20191201_v3.nc"

ds = fcpy.open_dataset(fpath)
ds = ds[["q"]]  # Select q only as this dataset contains more vars...
ds

<xarray.Dataset>
Dimensions:  (time: 1, lev: 137, lat: 451, lon: 900)
Coordinates:
  * time     (time) datetime64[ns] 2019-12-01T12:00:00
  * lon      (lon) float64 0.0 0.4 0.8 1.2 1.6 ... 358.0 358.4 358.8 359.2 359.6
  * lat      (lat) float64 -90.0 -89.6 -89.2 -88.8 -88.4 ... 88.8 89.2 89.6 90.0
  * lev      (lev) float64 1.0 2.0 3.0 4.0 5.0 ... 133.0 134.0 135.0 136.0 137.0
Data variables:
    q        (time, lev, lat, lon) float32 dask.array<chunksize=(1, 137, 451, 900), meta=np.ndarray>
Attributes:
    CDI:          Climate Data Interface version 1.9.8 (https://mpimet.mpg.de...
    Conventions:  CF-1.6
    history:      Mon Feb 15 18:37:28 2021: cdo -f nc4 copy tmp.grib data/mil...
    institution:  European Centre for Medium-Range Weather Forecasts
    CDO:          Climate Data Operators version 1.9.8 (https://mpimet.mpg.de...

xarray.Dataset

• Dimensions:
  • time: 1
  • lev: 137
  • lat: 451
  • lon: 900
• Coordinates: (4)
  • time
    (time)
    datetime64[ns]
    2019-12-01T12:00:00

    standard_name :

    time

    axis :

    T

    array(['2019-12-01T12:00:00.000000000'], dtype='datetime64[ns]')

  • lon
    (lon)
    float64
    0.0 0.4 0.8 ... 358.8 359.2 359.6

    standard_name :

    longitude

    long_name :

    longitude

    units :

    degrees_east

    axis :

    X

    array([  0. ,   0.4,   0.8, ..., 358.8, 359.2, 359.6])

  • lat
    (lat)
    float64
    -90.0 -89.6 -89.2 ... 89.6 90.0

    standard_name :

    latitude

    long_name :

    latitude

    units :

    degrees_north

    axis :

    Y

    array([-90. , -89.6, -89.2, ...,  89.2,  89.6,  90. ])

  • lev
    (lev)
    float64
    1.0 2.0 3.0 ... 135.0 136.0 137.0

    standard_name :

    hybrid_sigma_pressure

    long_name :

    hybrid level at layer midpoints

    formula :

    hyam hybm (mlev=hyam+hybm*aps)

    formula_terms :

    ap: hyam b: hybm ps: aps

    units :

    level

    positive :

    down

    array([  1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,   9.,  10.,  11.,  12.,
            13.,  14.,  15.,  16.,  17.,  18.,  19.,  20.,  21.,  22.,  23.,  24.,
            25.,  26.,  27.,  28.,  29.,  30.,  31.,  32.,  33.,  34.,  35.,  36.,
            37.,  38.,  39.,  40.,  41.,  42.,  43.,  44.,  45.,  46.,  47.,  48.,
            49.,  50.,  51.,  52.,  53.,  54.,  55.,  56.,  57.,  58.,  59.,  60.,
            61.,  62.,  63.,  64.,  65.,  66.,  67.,  68.,  69.,  70.,  71.,  72.,
            73.,  74.,  75.,  76.,  77.,  78.,  79.,  80.,  81.,  82.,  83.,  84.,
            85.,  86.,  87.,  88.,  89.,  90.,  91.,  92.,  93.,  94.,  95.,  96.,
            97.,  98.,  99., 100., 101., 102., 103., 104., 105., 106., 107., 108.,
           109., 110., 111., 112., 113., 114., 115., 116., 117., 118., 119., 120.,
           121., 122., 123., 124., 125., 126., 127., 128., 129., 130., 131., 132.,
           133., 134., 135., 136., 137.])

• Data variables: (1)
  • q
    (time, lev, lat, lon)
    float32
    dask.array<chunksize=(1, 137, 451, 900), meta=np.ndarray>

    standard_name :

    specific_humidity

    long_name :

    Specific humidity

    units :

    kg kg**-1

    param :

    0.1.0

    institution :

    ECMWF

    Array Chunk Bytes 212.13 MiB 212.13 MiB Shape (1, 137, 451, 900) (1, 137, 451, 900) Dask graph 1 chunks in 2 graph layers Data type float32 numpy.ndarray

• Indexes: (4)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2019-12-01 12:00:00'], dtype='datetime64[ns]', name='time', freq=None))

  • lon
    PandasIndex

    PandasIndex(Index([               0.0,                0.4,                0.8,
           1.2000000000000002,                1.6,                2.0,
           2.4000000000000004, 2.8000000000000003,                3.2,
                          3.6,
           ...
                        356.0, 356.40000000000003,              356.8,
           357.20000000000005,              357.6,              358.0,
           358.40000000000003,              358.8, 359.20000000000005,
                        359.6],
          dtype='float64', name='lon', length=900))

  • lat
    PandasIndex

    PandasIndex(Index([            -90.0,             -89.6,             -89.2,
                       -88.8,             -88.4,             -88.0,
                       -87.6,             -87.2,             -86.8,
                       -86.4,
           ...
                        86.4, 86.80000000000001, 87.20000000000002,
           87.60000000000002,              88.0,              88.4,
           88.80000000000001, 89.20000000000002, 89.60000000000002,
                        90.0],
          dtype='float64', name='lat', length=451))

  • lev
    PandasIndex

    PandasIndex(Index([  1.0,   2.0,   3.0,   4.0,   5.0,   6.0,   7.0,   8.0,   9.0,  10.0,
           ...
           128.0, 129.0, 130.0, 131.0, 132.0, 133.0, 134.0, 135.0, 136.0, 137.0],
          dtype='float64', name='lev', length=137))

• Attributes: (5)

  CDI :

  Climate Data Interface version 1.9.8 (https://mpimet.mpg.de/cdi)

  Conventions :

  CF-1.6

  history :

  Mon Feb 15 18:37:28 2021: cdo -f nc4 copy tmp.grib data/milan/gg/cams_q_20191201.grib_v3.nc

  institution :

  European Centre for Medium-Range Weather Forecasts

  CDO :

  Climate Data Operators version 1.9.8 (https://mpimet.mpg.de/cdo)

Define own metric

# Advance use: define your own metric.
# Here we count the number of unique values per compressor type
# without relying on auto-chunking or precomputed histograms
# this gives full flexibility but with greater complexity
def unique_count(chunks, baseline, compressors, bits):
    from fcpy import run_compressor_single
    import xarray as xr
    from collections import defaultdict

    # Unique values of decompressed dataset
    # ... over all compressors and bits
    unique = defaultdict(set)
    for chunk in chunks:
        for compressor in compressors:
            for bits_ in bits:
                da_decompressed = run_compressor_single(chunk, compressor, bits_)
                unique[(compressor.name, bits_)] |= set(
                    da_decompressed.values.flatten()
                )

    counts = xr.DataArray(
        0,
        dims=["compressor", "bits"],
        coords={"compressor": [c.name for c in compressors], "bits": bits},
    )
    for compressor in compressors:
        for bits_ in bits:
            counts.loc[dict(compressor=compressor.name, bits=bits_)] = len(
                unique[(compressor.name, bits_)]
            )
    return counts

Define and run experiment

suite = fcpy.Suite(
    ds,
    baseline=fcpy.Float(bits=32),
    compressors=[fcpy.LinQuantization(), fcpy.Round()],
    metrics=[fcpy.Difference, fcpy.AbsoluteError],
    custom_metrics=[unique_count],
    bits=[12, 14, 16, 18],
    max_chunk_size_bytes=451 * 900 * 4,
    skip_histograms=True,
)

  Activating project at `~/work/field-compression/field-compression`
WARNING: method definition for == at /usr/share/miniconda/envs/fcpy/share/julia/packages/ChainRulesCore/ctmSK/src/tangent_types/tangent.jl:68 declares type variable T but does not use it.
WARNING: method definition for getindex at /usr/share/miniconda/envs/fcpy/share/julia/packages/ChainRulesCore/ctmSK/src/tangent_types/tangent.jl:120 declares type variable T but does not use it.
WARNING: method definition for getindex at /usr/share/miniconda/envs/fcpy/share/julia/packages/ChainRulesCore/ctmSK/src/tangent_types/tangent.jl:120 declares type variable P but does not use it.
WARNING: method definition for canonicalize at /usr/share/miniconda/envs/fcpy/share/julia/packages/ChainRulesCore/ctmSK/src/tangent_types/tangent.jl:240 declares type variable L but does not use it.
WARNING: method definition for canonicalize at /usr/share/miniconda/envs/fcpy/share/julia/packages/ChainRulesCore/ctmSK/src/tangent_types/tangent.jl:241 declares type variable L but does not use it.

Plot results

# As data are in xarray, custom plots and comparisons are easy!
# Verbosity here is to showcase full customization
import matplotlib.pyplot as plt

ds_unique_count = suite.custom_metrics[0]
ds_unique_count.q.plot.line(x="bits")
plt.title(f"{ds.q.long_name} in {ds.q.units}")
plt.xlabel("Bits")
plt.ylabel("Number of Unique Values");