Programmatic example with MARS/GRIB input

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

Imports

import os
from pathlib import Path
import json
import numpy as np
import fcpy

Load data

if (Path().home() / ".ecmwfapirc").exists():
    fpath = "../samples/cases/hplp_ml_q.json"
    request = json.load(open(fpath))["mars"]
    ds = fcpy.fetch_mars_dataset(request)
else:
    ds = fcpy.open_dataset("../data/mars.grib")

ds = ds.isel(hybrid=[0, 1])
ds

In  , overriding the default value (cache=True) with cache=False is not recommended.

<xarray.Dataset>
Dimensions:     (time: 1, step: 1, hybrid: 2, values: 654400)
Coordinates:
  * time        (time) datetime64[ns] 2020-07-21
  * step        (step) timedelta64[ns] 00:00:00
  * hybrid      (hybrid) float64 1.0 2.0
    latitude    (values) float64 ...
    longitude   (values) float64 ...
    valid_time  (time, step) datetime64[ns] ...
Dimensions without coordinates: values
Data variables:
    q           (time, step, hybrid, values) float64 ...
Attributes:
    GRIB_edition:            2
    GRIB_centre:             ecmf
    GRIB_centreDescription:  European Centre for Medium-Range Weather Forecasts
    GRIB_subCentre:          0
    Conventions:             CF-1.7
    institution:             European Centre for Medium-Range Weather Forecasts
    history:                 2023-05-08T18:01 GRIB to CDM+CF via cfgrib-0.9.1...
    path:                    ../data/mars.grib

xarray.Dataset

• Dimensions:
  • time: 1
  • step: 1
  • hybrid: 2
  • values: 654400
• Coordinates: (6)
  • time
    (time)
    datetime64[ns]
    2020-07-21

    long_name :

    initial time of forecast

    standard_name :

    forecast_reference_time

    array(['2020-07-21T00:00:00.000000000'], dtype='datetime64[ns]')

  • step
    (step)
    timedelta64[ns]
    00:00:00

    long_name :

    time since forecast_reference_time

    standard_name :

    forecast_period

    array([0], dtype='timedelta64[ns]')

  • hybrid
    (hybrid)
    float64
    1.0 2.0

    long_name :

    hybrid level

    units :

    1

    positive :

    down

    standard_name :

    atmosphere_hybrid_sigma_pressure_coordinate

    array([1., 2.])

  • latitude
    (values)
    float64
    ...

    units :

    degrees_north

    standard_name :

    latitude

    long_name :

    latitude

    [654400 values with dtype=float64]

  • longitude
    (values)
    float64
    ...

    units :

    degrees_east

    standard_name :

    longitude

    long_name :

    longitude

    [654400 values with dtype=float64]

  • valid_time
    (time, step)
    datetime64[ns]
    ...

    standard_name :

    time

    long_name :

    time

    [1 values with dtype=datetime64[ns]]

• Data variables: (1)
  • q
    (time, step, hybrid, values)
    float64
    ...

    GRIB_paramId :

    133

    GRIB_dataType :

    fc

    GRIB_numberOfPoints :

    654400

    GRIB_typeOfLevel :

    hybrid

    GRIB_stepUnits :

    1

    GRIB_stepType :

    instant

    GRIB_gridType :

    reduced_gg

    GRIB_N :

    400

    GRIB_NV :

    276

    GRIB_cfName :

    specific_humidity

    GRIB_cfVarName :

    q

    GRIB_gridDefinitionDescription :

    Gaussian latitude/longitude

    GRIB_missingValue :

    9999

    GRIB_name :

    Specific humidity

    GRIB_pl :

    [ 20 24 28 32 36 40 44 48 52 56 60 64 68 72 76 80 84 88 92 96 100 104 108 112 116 120 124 128 132 136 140 144 148 152 156 160 164 168 172 176 180 184 188 192 196 200 204 208 212 216 220 224 228 232 236 240 244 248 252 256 260 264 268 272 276 280 284 288 292 296 300 304 308 312 316 320 324 328 332 336 340 344 348 352 356 360 364 368 372 376 380 384 388 392 396 400 404 408 412 416 420 424 428 432 436 440 444 448 452 456 460 464 468 472 476 480 484 488 492 496 500 504 508 512 516 520 524 528 532 536 540 544 548 552 556 560 564 568 572 576 580 584 588 592 596 600 604 608 612 616 620 624 628 632 636 640 644 648 652 656 660 664 668 672 676 680 684 688 692 696 700 704 708 712 716 720 724 728 732 736 740 744 748 752 756 760 764 768 772 776 780 784 788 792 796 800 804 808 812 816 820 824 828 832 836 840 844 848 852 856 860 864 868 872 876 880 884 888 892 896 900 904 908 912 916 920 924 928 932 936 940 944 948 952 956 960 964 968 972 976 980 984 988 992 996 1000 1004 1008 1012 1016 1020 1024 1028 1032 1036 1040 1044 1048 1052 1056 1060 1064 1068 1072 1076 1080 1084 1088 1092 1096 1100 1104 1108 1112 1116 1120 1124 1128 1132 1136 1140 1144 1148 1152 1156 1160 1164 1168 1172 1176 1180 1184 1188 1192 1196 1200 1204 1208 1212 1216 1220 1224 1228 1232 1236 1240 1244 1248 1252 1256 1260 1264 1268 1272 1276 1280 1284 1288 1292 1296 1300 1304 1308 1312 1316 1320 1324 1328 1332 1336 1340 1344 1348 1352 1356 1360 1364 1368 1372 1376 1380 1384 1388 1392 1396 1400 1404 1408 1412 1416 1420 1424 1428 1432 1436 1440 1444 1448 1452 1456 1460 1464 1468 1472 1476 1480 1484 1488 1492 1496 1500 1504 1508 1512 1516 1520 1524 1528 1532 1536 1540 1544 1548 1552 1556 1560 1564 1568 1572 1576 1580 1584 1588 1592 1596 1600 1604 1608 1612 1616 1616 1612 1608 1604 1600 1596 1592 1588 1584 1580 1576 1572 1568 1564 1560 1556 1552 1548 1544 1540 1536 1532 1528 1524 1520 1516 1512 1508 1504 1500 1496 1492 1488 1484 1480 1476 1472 1468 1464 1460 1456 1452 1448 1444 1440 1436 1432 1428 1424 1420 1416 1412 1408 1404 1400 1396 1392 1388 1384 1380 1376 1372 1368 1364 1360 1356 1352 1348 1344 1340 1336 1332 1328 1324 1320 1316 1312 1308 1304 1300 1296 1292 1288 1284 1280 1276 1272 1268 1264 1260 1256 1252 1248 1244 1240 1236 1232 1228 1224 1220 1216 1212 1208 1204 1200 1196 1192 1188 1184 1180 1176 1172 1168 1164 1160 1156 1152 1148 1144 1140 1136 1132 1128 1124 1120 1116 1112 1108 1104 1100 1096 1092 1088 1084 1080 1076 1072 1068 1064 1060 1056 1052 1048 1044 1040 1036 1032 1028 1024 1020 1016 1012 1008 1004 1000 996 992 988 984 980 976 972 968 964 960 956 952 948 944 940 936 932 928 924 920 916 912 908 904 900 896 892 888 884 880 876 872 868 864 860 856 852 848 844 840 836 832 828 824 820 816 812 808 804 800 796 792 788 784 780 776 772 768 764 760 756 752 748 744 740 736 732 728 724 720 716 712 708 704 700 696 692 688 684 680 676 672 668 664 660 656 652 648 644 640 636 632 628 624 620 616 612 608 604 600 596 592 588 584 580 576 572 568 564 560 556 552 548 544 540 536 532 528 524 520 516 512 508 504 500 496 492 488 484 480 476 472 468 464 460 456 452 448 444 440 436 432 428 424 420 416 412 408 404 400 396 392 388 384 380 376 372 368 364 360 356 352 348 344 340 336 332 328 324 320 316 312 308 304 300 296 292 288 284 280 276 272 268 264 260 256 252 248 244 240 236 232 228 224 220 216 212 208 204 200 196 192 188 184 180 176 172 168 164 160 156 152 148 144 140 136 132 128 124 120 116 112 108 104 100 96 92 88 84 80 76 72 68 64 60 56 52 48 44 40 36 32 28 24 20]

    GRIB_shortName :

    q

    GRIB_units :

    kg kg**-1

    long_name :

    Specific humidity

    units :

    kg kg**-1

    standard_name :

    specific_humidity

    [1308800 values with dtype=float64]

• Indexes: (3)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2020-07-21'], dtype='datetime64[ns]', name='time', freq=None))

  • step
    PandasIndex

    PandasIndex(TimedeltaIndex(['0 days'], dtype='timedelta64[ns]', name='step', freq=None))

  • hybrid
    PandasIndex

    PandasIndex(Index([1.0, 2.0], dtype='float64', name='hybrid'))

• Attributes: (8)

  GRIB_edition :

  2

  GRIB_centre :

  ecmf

  GRIB_centreDescription :

  European Centre for Medium-Range Weather Forecasts

  GRIB_subCentre :

  0

  Conventions :

  CF-1.7

  institution :

  European Centre for Medium-Range Weather Forecasts

  history :

  2023-05-08T18:01 GRIB to CDM+CF via cfgrib-0.9.10.1/ecCodes-2.30.0 with {"source": "N/A", "filter_by_keys": {}, "encode_cf": ["parameter", "time", "geography", "vertical"]}

  path :

  ../data/mars.grib

Define and run experiment

suite = fcpy.Suite(
    ds,
    baseline=fcpy.Float(bits=32),
    compressors=[fcpy.LinQuantization(), fcpy.Round()],
    metrics=[fcpy.Difference, fcpy.AbsoluteError],
    bits=[13, 14, 15, 16, 17, 18, 19, 20, 21, 22],
    max_chunk_size_bytes=451 * 900 * 6,
    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.

Lineplot: evaluate the effect of bits

suite.lineplot(fcpy.AbsoluteError, reduction="max")

Boxplot: evaluate statistical characteristics of errors

bits = 13
suite.boxplot(metric=fcpy.Difference, compressor="LinQuantization", bits=bits)
suite.boxplot(metric=fcpy.Difference, compressor="Round", bits=bits)

Spatial plots: evaluate the spatial distribution of errors

if not os.getenv("SPHINX_BUILD"):
    # Here we compare the spatial change of specific humidity
    # by rounding values to 8 bits
    fcpy.spatialplot(
        ds=ds,
        baseline=fcpy.Float(bits=32),
        var_name="q",
        compressor=fcpy.Round(bits=6),
        metric=fcpy.Difference,
        latitude=0,
        longitude=0,
        hybrid=1,
        step=np.timedelta64(0, "ns"),
        third_dim="hybrid",
    )