diff --git a/.buildinfo b/.buildinfo index c1b61c9..dee73c8 100644 --- a/.buildinfo +++ b/.buildinfo @@ -1,4 +1,4 @@ # Sphinx build info version 1 # This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done. -config: e55ff71d5e34c542c153cb9972ccded6 +config: 7e2da3bb54b899b03b4b4974f7b4ab4b tags: 645f666f9bcd5a90fca523b33c5a78b7 diff --git a/CONTRIBUTING.html b/CONTRIBUTING.html index d55fac9..c408a3e 100644 --- a/CONTRIBUTING.html +++ b/CONTRIBUTING.html @@ -9,7 +9,7 @@ - Developer Guide for the website — ClimSim + Contributor Guidelines — ClimSim @@ -64,6 +64,7 @@ + @@ -98,8 +99,8 @@ class="form-control" name="q" id="search-input" - placeholder="Search..." - aria-label="Search..." + placeholder="Search this book..." + aria-label="Search this book..." autocomplete="off" autocorrect="off" autocapitalize="off" @@ -129,7 +130,7 @@
-
@@ -320,624 +329,45 @@

Contents

Dataset: E3SM-MMF Low-Resolution Aquaplanet

ClimSim: An open large-scale dataset for training high-resolution physics emulators in hybrid multi-scale climate simulators#

-

This repository contains the code necessary to download and preprocess the data, and create, train, and evaluate the baseline models in the paper.

+

Code Repository

+

ClimSim is the largest-ever dataset designed for hybrid ML-physics research. It comprises multi-scale climate simulations, developed by a consortium of climate scientists and ML researchers. It consists of 5.7 billion pairs of multivariate input and output vectors that isolate the influence of locally-nested, high-resolution, high-fidelity physics on a host climate simulator’s macro-scale physical state. +The dataset is global in coverage, spans multiple years at high sampling frequency, and is designed such that resulting emulators are compatible with downstream coupling into operational climate simulators. +We implement a range of deterministic and stochastic regression baselines to highlight the ML challenges and their scoring.

fig_1

-
-

Introductory video:#

-

https://www.youtube.com/watch?v=M3Vz0zR1Auc

-

Alt text

+
+

Project structure#

-
-

Dataset Information#

-

Data from multi-scale climate model (E3SM-MMF) simulations were saved at 20-minute intervals for 10 simulated years. Two netCDF files–input and output (target)–are produced at each timestep, totaling 525,600 files for each configuration. 3 configurations of E3SM-MMF were run:

-
    -

  1. High-Resolution Real Geography

    -
      -

    • 1.5° x 1.5° horizontal resolution (21,600 grid columns)

      • -

    • 5.7 billion total samples (41.2 TB)

      • -

    • 102 MB per input file, 61 MB per output file

      • -
    -
    2. -

  2. Low-Resolution Real Geography

    +
    +

    Getting Started#

    -
    3. -

  3. Low-Resolution Aquaplanet

    +
+
+

Demo notebooks#

- - -

Dataset used for baseline models showcased in paper corresponds to the Low-Resolution Real Geography dataset. Scalar variables vary in time and “horizontal” grid (“ncol”), while vertically-resolved variables vary additionally in vertical space (“lev”). The full list of variables can be found in Supplmentary Information Table 1. The subset of variables used in our experiments is shown below:

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

Input

Target

Variable

Description

Units

Dimensions

X

T

Air temperature

K

(lev, ncol)

X

q

Specific humidity

kg/kg

(lev, ncol)

X

PS

Surface pressure

Pa

(ncol)

X

SOLIN

Solar insolation

W/m²

(ncol)

X

LHFLX

Surface latent heat flux

W/m²

(ncol)

X

SHFLX

Surface sensible heat flux

W/m²

(ncol)

X

dT/dt

Heating tendency

K/s

(lev, ncol)

X

dq/dt

Moistening tendency

kg/kg/s

(lev, ncol)

X

NETSW

Net surface shortwave flux

W/m²

(ncol)

X

FLWDS

Downward surface longwave flux

W/m²

(ncol)

X

PRECSC

Snow rate

m/s

(ncol)

X

PRECC

Rain rate

m/s

(ncol)

X

SOLS

Visible direct solar flux

W/m²

(ncol)

X

SOLL

Near-IR direct solar flux

W/m²

(ncol)

X

SOLSD

Visible diffuse solar flux

W/m²

(ncol)

X

SOLLD

Near-IR diffuse solar flux

W/m²

(ncol)

-
-

Quickstart#

-

This method is intended for those who wish to quickly try out new methods on a manageable subset of the data uploaded to HuggingFace.

-

Step 1

-

The first step is to download the subsampled low-resolution real-geography version of the data here. This contains subsampled and prenormalized data that was used for training, validation, and metrics for the ClimSim paper. It can be reproduced with the full version of the dataset using the preprocessing/create_npy_data_splits.ipynb notebook.

-

Training data corresponds to train_input.npy and train_target.npy. Validation data corresponds to val_input.npy and val_target.npy. Scoring data (which can be treated as a test set) corresponds to scoring_input.npy and scoring_target.npy. We have an additional held-out test set that we will use for an upcoming online competition. Keep an eye out! 😉

-

Step 2

-

Install the climsim_utils python tools, by running the following code from the root of this repo:

-
pip install .
-
-
-

If you already have all climsim_utils dependencies (tensorflow, xarray, etc.) installed in your local environment, you can alternatively run:

-
pip install . --no-deps
-
-
-

Step 3

-

Train your model on the training data and validate using the validation data. If you wish to use something like a CNN, you will probably want to separate the variables into channels and broadcast scalars into vectors of the same dimension as vertically-resolved variables. Methods to do this can be found in the climsim_utils/data_utils.py script.

-

Step 4

-

Evaluation time! Use the evaluation/main_figure_generation.ipynb notebook to see how your model does! Use the calc_MAE, calc_RMSE, and calc_R2 methods in the climsim_utils/data_utils.py script to see how your model does on point estimates and use the calc_CRPS method to check how well-calibrated your model is if it’s stochastic. 😊

-
-
-

Download the Data#

-

The input (“mli”) and target (“mlo”) data for all E3SM-MMF configurations can be downloaded from Hugging Face:

+
+

References#

-

There is also a quickstart dataset that contains subsampled and prenormalized data. This data was used for training, validation, and metrics for the ClimSim paper and can be reproduced using the preprocessing/create_npy_data_splits.ipynb notebook.

+

Alt text

-
-

Installation & setup#

-

For preprocessing and evaluation, please install the climsim_utils python tools, by running the following code from the root of this repo:

-
pip install .
-
-
-

If you already have all climsim_utils dependencies (tensorflow, xarray, etc.) installed in your local environment, you can alternatively run:

-
pip install . --no-deps
-
+
-
-
-

Preprocess the Data#

-

The default preprocessing workflow takes folders of monthly data from the climate model simulations, and creates normalized NumPy arrays for input and target data for training, validation, and scoring. These NumPy arrays are called train_input.npy, train_target.npy, val_input.npy, val_target.npy, scoring_input.npy, and scoring_target.npy. An option to strictly use a data loader and avoid converting into NumPy arrays is available in data_utils.py; however, this can slow down training because of increased I/O.

-

The data comes in the form of folders labeled YYYY-MM, which corresponds to the simulation year (YYYY) and month (MM). Within each of these folders are netCDF (.nc) files that represent inputs and targets for individual timesteps. Input files are labeled E3SM-MMF.mli.YYYY-MM-DD-SSSSS.nc where DD-SSSSS corresponds to the day of the month (DD) and seconds of the day (SSSSS), with timesteps being spaced 1,200 seconds (20 minutes) apart. Target files are labeled the same way, except mli is replaced by mlo. For vertically-resolved variables, lower indices corresponds to higher levels in the atmosphere. This is because pressure decreases monotonically with altitude.

-

The files containing the default normalization factors for the input and target data are found in the normalizations/ folder, precomputed for convenience. However, one can use their own normalization factors if desired. The file containing the E3SM-MMF grid information is found in the grid_info/ folder. This corresponds to the netCDF file ending in grid-info.nc on Hugging Face.

-

The environment needed for preprocessing can be found in the /preprocessing/env/requirements.txt file. A class designed for preprocessing and metrics can be imported from the data_utils.py script. This script is used in the preprocessing/create_npy_data_splits.ipynb notebook, which creates training, validation, and scoring datasets.

-

By default, training and validation data subsample every \(7^{th}\) timestep while scoring data subsamples every \(6^{th}\) timestep to enable daily-averaged metrics. Training data is taken from the second month of simulation year 1 through the first month of simulation year 8 (i.e., 0001-02 through 0008-01). Both validation and scoring data are taken from 0008-02 through 0009-01. However, the data_utils.py allows the user to easily change these defaults assuming knowledge of regular expressions. To see how this works, please reference preprocessing/create_npy_data_splits.ipynb.

-
-
-

Baseline Models#

-

Six different baseline models were created and trained:

-
    -

  1. Convolutional neural network (CNN)

    2. -

  2. Encoder-decoder (ED)

    3. -

  3. Heteroskedastic regression (HSR)

    4. -

  4. Multi-layer perceptron (MLP)

    5. -

  5. Randomized prior network (RPN)

    6. -

  6. Conditional variational autoencoder (cVAE)

    7. -
-

Jupyter Notebooks describing how to load and train simple CNN and MLP models are found in the demo_notebooks/ folder. The environments and code used to train each model, as well as the pre-trained models, are found in the baseline_models/ folder.

-
-
-

Evaluation#

-

Four different evaluation metrics were calculated:

-
    -

  1. Mean absolute error (MAE)

    2. -

  2. Coefficient of determination (R²)

    3. -

  3. Root mean squared error (RMSE)

    4. -

  4. Continuous ranked probability score (CRPS)

    5. -
-

Evaluation and comparison of the different baseline models are found in the evaluation/ folder. All variables are converted to a common energy unit (i.e., W/m²) for scoring. The scoring is done using the functions in evaluation/data_utils.py.

-

Evaluation metrics are computed separately for each horizontally-averaged, vertically-averaged, and time-averaged target variable. The performance for each baseline model for all four metrics is shown below:

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

MAE (W/m²)

CNN

ED

HSR

MLP

RPN

cVAE

dT/dt

2.585

2.684

2.845

2.683

2.685

2.732

dq/dt

4.401

4.673

4.784

4.495

4.592

4.680

NETSW

18.85

14.968

19.82

13.36

18.88

19.73

FLWDS

8.598

6.894

6.267

5.224

6.018

6.588

PRECSC

3.364

3.046

3.511

2.684

3.328

3.322

PRECC

37.83

37.250

42.38

34.33

37.46

38.81

SOLS

10.83

8.554

11.31

7.97

10.36

10.94

SOLL

13.15

10.924

13.60

10.30

12.96

13.46

SOLSD

5.817

5.075

6.331

4.533

5.846

6.159

SOLLD

5.679

5.136

6.215

4.806

5.702

6.066

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

R²

CNN

ED

HSR

MLP

RPN

cVAE

dT/dt

0.627

0.542

0.568

0.589

0.617

0.590

dq/dt

–

–

–

–

–

–

NETSW

0.944

0.980

0.959

0.983

0.968

0.957

FLWDS

0.828

0.802

0.904

0.924

0.912

0.883

PRECSC

–

–

–

–

–

–

PRECC

0.077

-17.909

-68.35

-38.69

-67.94

-0.926

SOLS

0.927

0.960

0.929

0.961

0.943

0.929

SOLL

0.916

0.945

0.916

0.948

0.928

0.915

SOLSD

0.927

0.951

0.923

0.956

0.940

0.921

SOLLD

0.813

0.857

0.797

0.866

0.837

0.796

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

RMSE (W/m²)

CNN

ED

HSR

MLP

RPN

cVAE

dT/dt

4.369

4.696

4.825

4.421

4.482

4.721

dq/dt

7.284

7.643

7.896

7.322

7.518

7.780

NETSW

36.91

28.537

37.77

26.71

33.60

38.36

FLWDS

10.86

9.070

8.220

6.969

7.914

8.530

PRECSC

6.001

5.078

6.095

4.734

5.511

6.182

PRECC

85.31

76.682

90.64

72.88

76.58

88.71

SOLS

22.92

17.999

23.61

17.40

20.61

23.27

SOLL

27.25

22.540

27.78

21.95

25.22

27.81

SOLSD

12.13

9.917

12.40

9.420

11.00

12.64

SOLLD

12.10

10.417

12.47

10.12

11.25

12.63

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

CRPS (W/m²)

CNN

ED

HSR

MLP

RPN

cVAE

dT/dt

–

–

3.284

–

2.580

2.795

dq/dt

–

–

4.899

–

4.022

4.372

NETSW

–

–

0.055

–

0.053

0.057

FLWDS

–

–

0.018

–

0.016

0.018

PRECSC

–

–

0.011

–

0.008

0.009

PRECC

–

–

0.122

–

0.085

0.097

SOLS

–

–

0.031

–

0.028

0.033

SOLL

–

–

0.038

–

0.035

0.040

SOLSD

–

–

0.018

–

0.015

0.016

SOLLD

–

–

0.017

–

0.015

0.016

-

The evaluation/main_figure_generation.ipynb notebook calculates and plots MAE, R², RMSE, and CRPS scores for each baseline model. The separate R² for longitudinally-averaged and time-averaged 3D variables is found in plot_R2_analysis.ipynb.

-
-
-

Feedback and future development#

-

Additional issues from community researchers wishing to build off ClimSim or make use of our tools are welcome and can be raised using the GitHub issues page (preferred) or by directly emailing ClimSim maintainers.

-
- + @@ -100,8 +100,8 @@ class="form-control" name="q" id="search-input" - placeholder="Search..." - aria-label="Search..." + placeholder="Search this book..." + aria-label="Search this book..." autocomplete="off" autocorrect="off" autocapitalize="off" @@ -131,7 +131,7 @@
-
+ @@ -101,8 +102,8 @@ class="form-control" name="q" id="search-input" - placeholder="Search..." - aria-label="Search..." + placeholder="Search this book..." + aria-label="Search this book..." autocomplete="off" autocorrect="off" autocapitalize="off" @@ -132,7 +133,7 @@
-
+ +
+

next

+

Contributor Guidelines

+
+ + diff --git a/evaluating.html b/evaluating.html new file mode 100644 index 0000000..d725fe9 --- /dev/null +++ b/evaluating.html @@ -0,0 +1,814 @@ + + + + + + + + + + + + Evaluation — ClimSim + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+
+
+ + + Ctrl+K +
+
+ + + +
+
+ +
+ + + +
+ + + + +
+ +
+
+ + +
+
+
+ + +
+
+ + + + +
+ +
+ + + +
+ +
+
+ +
+
+ +
+ +
+ +
+ + +
+ +
+ +
+ + + + + +
+ + +
+ + + + + + + + + + + +
+ +
+ +
+
+ + + +
+

Evaluation

+ +
+
+ +
+
+
+ + + + +
+ +
+

Evaluation#

+

Four different evaluation metrics were calculated:

+
    +

  1. Mean absolute error (MAE)

    2. +

  2. Coefficient of determination (R²)

    3. +

  3. Root mean squared error (RMSE)

    4. +

  4. Continuous ranked probability score (CRPS)

    5. +
+

Evaluation and comparison of the different baseline models are found in the [evaluation/](https://github.com/leap-stc/ClimSim/tree/main/evaluation) folder on GitHub. All variables are converted to a common energy unit (i.e., W/m²) for scoring. The scoring is done using the functions in [climsim_utils/data_utils.py](https://github.com/leap-stc/ClimSim/tree/main/climsim_utils).

+

This notebook calculates and plots MAE, R², RMSE, and CRPS scores for each baseline model. The separate R² for longitudinally-averaged and time-averaged 3D variables is found in this notebook.

+

Evaluation metrics are computed separately for each horizontally-averaged, vertically-averaged, and time-averaged target variable. The performance for each baseline model for all four metrics is shown below:

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +

MAE (W/m²)

CNN

ED

HSR

MLP

RPN

cVAE

dT/dt

2.585

2.684

2.845

2.683

2.685

2.732

dq/dt

4.401

4.673

4.784

4.495

4.592

4.680

NETSW

18.85

14.968

19.82

13.36

18.88

19.73

FLWDS

8.598

6.894

6.267

5.224

6.018

6.588

PRECSC

3.364

3.046

3.511

2.684

3.328

3.322

PRECC

37.83

37.250

42.38

34.33

37.46

38.81

SOLS

10.83

8.554

11.31

7.97

10.36

10.94

SOLL

13.15

10.924

13.60

10.30

12.96

13.46

SOLSD

5.817

5.075

6.331

4.533

5.846

6.159

SOLLD

5.679

5.136

6.215

4.806

5.702

6.066

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +

R²

CNN

ED

HSR

MLP

RPN

cVAE

dT/dt

0.627

0.542

0.568

0.589

0.617

0.590

dq/dt

–

–

–

–

–

–

NETSW

0.944

0.980

0.959

0.983

0.968

0.957

FLWDS

0.828

0.802

0.904

0.924

0.912

0.883

PRECSC

–

–

–

–

–

–

PRECC

0.077

-17.909

-68.35

-38.69

-67.94

-0.926

SOLS

0.927

0.960

0.929

0.961

0.943

0.929

SOLL

0.916

0.945

0.916

0.948

0.928

0.915

SOLSD

0.927

0.951

0.923

0.956

0.940

0.921

SOLLD

0.813

0.857

0.797

0.866

0.837

0.796

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +

RMSE (W/m²)

CNN

ED

HSR

MLP

RPN

cVAE

dT/dt

4.369

4.696

4.825

4.421

4.482

4.721

dq/dt

7.284

7.643

7.896

7.322

7.518

7.780

NETSW

36.91

28.537

37.77

26.71

33.60

38.36

FLWDS

10.86

9.070

8.220

6.969

7.914

8.530

PRECSC

6.001

5.078

6.095

4.734

5.511

6.182

PRECC

85.31

76.682

90.64

72.88

76.58

88.71

SOLS

22.92

17.999

23.61

17.40

20.61

23.27

SOLL

27.25

22.540

27.78

21.95

25.22

27.81

SOLSD

12.13

9.917

12.40

9.420

11.00

12.64

SOLLD

12.10

10.417

12.47

10.12

11.25

12.63

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +

CRPS (W/m²)

CNN

ED

HSR

MLP

RPN

cVAE

dT/dt

–

–

3.284

–

2.580

2.795

dq/dt

–

–

4.899

–

4.022

4.372

NETSW

–

–

0.055

–

0.053

0.057

FLWDS

–

–

0.018

–

0.016

0.018

PRECSC

–

–

0.011

–

0.008

0.009

PRECC

–

–

0.122

–

0.085

0.097

SOLS

–

–

0.031

–

0.028

0.033

SOLL

–

–

0.038

–

0.035

0.040

SOLSD

–

–

0.018

–

0.015

0.016

SOLLD

–

–

0.017

–

0.015

0.016

+
+
+
+ + + + +
+ + + + + + +
+ + + + +
+
+ +
+ +
+ +

+By The ClimSim Maintainers +

+ +
+ +
+ +

+ + © Copyright 2023. +
+ +

+ +
+ +
+ +
+ +
+ +
+ +
+
+ + +
+
+
+ + + + + + + + \ No newline at end of file diff --git a/evaluation/ClimSim_metrics.html b/evaluation/ClimSim_metrics.html new file mode 100644 index 0000000..8681e39 --- /dev/null +++ b/evaluation/ClimSim_metrics.html @@ -0,0 +1,1887 @@ + + + + + + + + + + + + ClimSim: Metrics calculation and visualization — ClimSim + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+
+
+ + + Ctrl+K +
+
+ + + +
+
+ +
+ + + +
+ + + + +
+ +
+
+ + +
+
+
+ + +
+
+ + + + +
+ +
+ + + +
+ +
+
+ +
+
+ +
+ +
+ +
+ + +
+ +
+ +
+ + + + + +
+ + +
+ + + + + + + + + + + +
+ +
+ +
+
+ + + +
+

ClimSim: Metrics calculation and visualization

+ +
+ +
+
+ + + + +
+ +
+

ClimSim: Metrics calculation and visualization#

+

Author: Sungduk Yu
+Last update: Wed 16 Aug 2023 04:52:19 PM PDT
+Description:

+
    +

  • This script calculates evaluation three metrics (MAE, R2, RMSE) in Part 1 and generates tables and figures for ClimSim manuscript in Part 2.\

    • +

  • CRPS is calculated by a separate python script (“crps_clean.py” written by bhouri0412). In Part 2, presaved CRPS metrics as npy files are read and visualized.

    • +
+

Input:

+
    +

  • model_name: a string, model name

    • +

  • fn_x_true: npy file name for true input

    • +

  • fn_y_true: npy file name for true output

    • +

  • fn_y_pred: npy file name for model prediction

    • +

  • fn_grid: netcdif file name for grid information

    • +

  • fn_mli_mean, fn_mli_min, fn_mli_max: netcdf file names for input normalization coefficients (mean, min, max)

    • +

  • fn_mlo_scale: netcdf file name for output scaling coefficients.

    • +

  • fn_save_metrics: csv file name for saving output

    • +

  • fn_save_metrics_avg: csv file name for saving output (avg over the level dimension)

    • +
+

Output:

+
    +

  • calculated metrics as csv file with names set in fn_save_metrics and fn_save_metrics_avg

    • +
+
+
+

<PART 1> Calculate metrics and save them into a csv file#

+
+

Change file names [user input here only!]#

+
+
+
FN_MODEL_OUTPUT = {'MLP':  './model_outputs/001_backup_phase-7_retrained_models_step2_lot-147_trial_0027.best.h5.npy',
+                   'RPN':  './model_outputs/rpn_pred_v1_stride6.npy',
+                   'CNN':  './model_outputs/val_predict_cnn_reshaped_stride6_FINAL.npy',
+                   'cVAE': './model_outputs/cvae.h5',
+                   'HSR': './model_outputs/hsr_preds_bestcrps.h5',
+                   'ED': './model_outputs/ED_ClimSIM_1_3_pred.npy'
+                  }
+
+
+
+
+
+
+
# model name
+# (model name is used for the output)
+model_name = 'MLP'
+
+# input of validation dataset (npy)
+fn_x_true = '../npy_data_conversion/npy_files/val_input_stride6.npy'
+
+# true output of validation dataset (npy)
+fn_y_true = '../npy_data_conversion/npy_files/val_target_stride6.npy'
+
+# Model predicted output of varlidation dataset (npy)
+fn_y_pred = FN_MODEL_OUTPUT[model_name]
+
+# model grid information (nc)
+fn_grid = '../grid_info/E3SM-MMF_ne4_grid-info.orig.nc'
+
+# normalization scale factors (nc)
+fn_mli_mean  = '../norm_factors/mli_mean.nc'
+fn_mli_min   = '../norm_factors/mli_min.nc'
+fn_mli_max   = '../norm_factors/mli_max.nc'
+fn_mlo_scale = '../norm_factors/mlo_scale.nc'
+
+# fn_save_output
+fn_save_metrics = f'./metrics/{model_name}.metrics.csv'
+fn_save_metrics_avg = f'./metrics/{model_name}.metrics.lev-avg.csv'
+
+
+
+
+
+
+

Set energy conversion scale factors#

+

We want to convert all output variables to a common energy unit, e.g., W/m2.

+

Energy unit conversion

+
    +

  • “dT/dt” [K/s] * cp [J/kg/K] * dp/g [kg/m2] -> [W/m2]

    • +

  • “dQ/dt” [kg/kg/s] * lv [J/kg] * dp/g [kg/m2] -> [W/m2]

    • +

  • “PRECC” [m/s] * lv [J/kg] * rho_h2o [kg/m3] -> [W/m2]

    • +
+

(for dT/dt and dQ/dt, it should be “dp/g-weighted” vertical integration, not vertical averaging.)

+
+
+
# physical constatns from (E3SM_ROOT/share/util/shr_const_mod.F90)
+grav    = 9.80616    # acceleration of gravity ~ m/s^2
+cp      = 1.00464e3  # specific heat of dry air   ~ J/kg/K
+lv      = 2.501e6    # latent heat of evaporation ~ J/kg
+lf      = 3.337e5    # latent heat of fusion      ~ J/kg
+ls      = lv + lf    # latent heat of sublimation ~ J/kg
+rho_air = 101325./ (6.02214e26*1.38065e-23/28.966) / 273.15 # density of dry air at STP  ~ kg/m^3
+                                                            # ~ 1.2923182846924677
+                                                            # SHR_CONST_PSTD/(SHR_CONST_RDAIR*SHR_CONST_TKFRZ)
+                                                            # SHR_CONST_RDAIR   = SHR_CONST_RGAS/SHR_CONST_MWDAIR
+                                                            # SHR_CONST_RGAS    = SHR_CONST_AVOGAD*SHR_CONST_BOLTZ
+rho_h20 = 1.e3       # density of fresh water     ~ kg/m^ 3
+
+vars_mlo_energy_conv = {'ptend_t':cp,
+                        'ptend_q0001':lv,
+                        'cam_out_NETSW':1.,
+                        'cam_out_FLWDS':1.,
+                        'cam_out_PRECSC':lv*rho_h20,
+                        'cam_out_PRECC':lv*rho_h20,
+                        'cam_out_SOLS':1.,
+                        'cam_out_SOLL':1.,
+                        'cam_out_SOLSD':1.,
+                        'cam_out_SOLLD':1.
+                       }
+vars_longname=\
+{'ptend_t':'Heating tendency, ∂T/∂t',
+ 'ptend_q0001':'Moistening tendency, ∂q/∂t',
+ 'cam_out_NETSW':'Net surface shortwave flux, NETSW',
+ 'cam_out_FLWDS':'Downward surface longwave flux, FLWDS',
+ 'cam_out_PRECSC':'Snow rate, PRECSC',
+ 'cam_out_PRECC':'Rain rate, PRECC',
+ 'cam_out_SOLS':'Visible direct solar flux, SOLS',
+ 'cam_out_SOLL':'Near-IR direct solar flux, SOLL',
+ 'cam_out_SOLSD':'Visible diffused solar flux, SOLSD',
+ 'cam_out_SOLLD':'Near-IR diffused solar flux, SOLLD'}
+
+
+
+
+
+
+

Main#

+
+
+
import numpy as np
+import xarray as xr
+import pandas as pd
+import matplotlib.pyplot as plt
+import glob
+
+# set dimemsion names for xarray datasets
+dim_name_level  = 'lev'
+dim_name_sample = 'sample'
+
+
+
+
+
+
+
# load input dataset
+x_true = np.load(fn_x_true).astype(np.float64)
+y_true = np.load(fn_y_true).astype(np.float64)
+if fn_y_pred[-3:] == '.h5':
+    y_pred = xr.open_dataset(fn_y_pred)['pred'].values
+else:
+    y_pred = np.load(fn_y_pred).astype(np.float64)
+N_samples = y_pred.shape[0] 
+
+# load norm/scale factors
+mlo_scale = xr.open_dataset(fn_mlo_scale)
+mli_mean  = xr.open_dataset(fn_mli_mean)
+mli_min   = xr.open_dataset(fn_mli_min)
+mli_max   = xr.open_dataset(fn_mli_max)
+
+
+
+
+
+
+
# load grid information
+ds_grid = xr.open_dataset(fn_grid) # has ncol:384
+N_ncol = len(ds_grid['ncol']) # length of ncol dimension (nlat * nlon)
+
+# make area-weights
+ds_grid['area_wgt'] = ds_grid['area'] / ds_grid['area'].mean('ncol')
+
+# map ds_grid's ncol dimension -> the N_samples dimension of npy-loayd arrays (e.g., y_pred)
+to_xarray = {'area_wgt': (dim_name_sample,np.tile(ds_grid['area_wgt'], int(N_samples/len(ds_grid['ncol'])))),
+            }
+to_xarray = xr.Dataset(to_xarray)
+
+# add nsample-mapped grid variables back to ds_grid
+ds_grid = xr.merge([ds_grid  [['P0', 'hyai', 'hyam','hybi','hybm','lat','lon','area']],
+                    to_xarray[['area_wgt']]])
+
+
+
+
+
+

Pack np arrays to xarray dataset#

+

(This is only to use Xarray for metric calculations.)

+
+
+
# list of ML output variables
+vars_mlo = ['ptend_t','ptend_q0001','cam_out_NETSW','cam_out_FLWDS','cam_out_PRECSC',
+            'cam_out_PRECC','cam_out_SOLS','cam_out_SOLL','cam_out_SOLSD','cam_out_SOLLD'] # mlo mean ML output.
+
+# length of each variable
+# (make sure that the order of variables are correct)
+vars_mlo_len = {'ptend_t':60,
+                'ptend_q0001':60,
+                'cam_out_NETSW':1,
+                'cam_out_FLWDS':1,
+                'cam_out_PRECSC':1,
+                'cam_out_PRECC':1,
+                'cam_out_SOLS':1,
+                'cam_out_SOLL':1,
+                'cam_out_SOLSD':1,
+                'cam_out_SOLLD':1
+               }
+
+# map the length of dimension to the name of dimension
+len_to_dim = {60:dim_name_level,
+              N_samples: dim_name_sample}
+
+
+
+
+
+
+
# Here, we first construct a dictionary of {var name: (dimension name, array-like)},
+# then, map the dictionary to an Xarray Dataset.
+# (ref: https://docs.xarray.dev/en/stable/generated/xarray.Dataset.html)
+
+DS = {}
+
+for kds in ['true', 'pred']:
+    if kds=='true':
+        work = y_true
+    elif kds=='pred':
+        work = y_pred
+
+    # [1] Construct dictionary for xarray dataset
+    #     format is key for variable name /
+    #               value for a turple of (dimension names, data).
+    to_xarray = {}
+    for k, kvar in enumerate(vars_mlo):
+
+        # length of variable (ie, number of levels)
+        kvar_len = vars_mlo_len[kvar]
+
+        # set dimensions of variable
+        if kvar_len == 60:
+            kvar_dims = (dim_name_sample, dim_name_level)
+        elif kvar_len == 1:
+            kvar_dims = dim_name_sample
+
+        # set start and end indices of variable in the loaded numpy array
+        # then, add 'kvar':(kvar_dims, <np_array>) to dictionary
+        if k==0: ind1=0
+        ind2 = ind1 + kvar_len
+
+        # scaled output
+        kvar_data = np.squeeze(work[:,ind1:ind2])
+        # unscaled output
+        kvar_data = kvar_data / mlo_scale[kvar].values
+
+        to_xarray[kvar] = (kvar_dims, kvar_data)
+
+        ind1 = ind2
+
+    # [2] convert dict to xarray dataset
+    DS[kds] = xr.Dataset(to_xarray)
+
+    # [3] add surface pressure ('state_ps') from ml input
+    # normalized ps
+    state_ps =  xr.DataArray(x_true[:,120], dims=('sample'), name='state_ps')
+    # denormalized ps
+    state_ps = state_ps * (mli_max['state_ps'] - mli_min['state_ps']) + mli_mean['state_ps']
+    DS[kds]['state_ps'] = state_ps
+
+    # [4] add grid information
+    DS[kds] = xr.merge([DS[kds], ds_grid])
+
+    # [5] add pressure thickness of each level, dp
+    # FYI, in a hybrid sigma vertical coordinate system, pressure at level z is
+    # P[x,z] = hyam[z]*P0 + hybm[z]*PS[x,z],
+    # where, hyam and hybm are 
+    tmp = DS[kds]['P0']*DS[kds]['hyai'] + DS[kds]['state_ps']*DS[kds]['hybi']
+    tmp = tmp.isel(ilev=slice(1,61)).values - tmp.isel(ilev=slice(0,60)).values
+    tmp = tmp.transpose()
+    DS[kds]['dp'] = xr.DataArray(tmp, dims=('sample', 'lev'))
+
+    # [6] break (sample) to (ncol,time)
+    N_timestep = int(N_samples/N_ncol)
+    dim_ncol     = np.arange(N_ncol)
+    dim_timestep = np.arange(N_timestep)
+    new_ind = pd.MultiIndex.from_product([dim_timestep, dim_ncol],
+                                         names=['time', 'ncol'])
+    DS[kds] = DS[kds].assign_coords(sample=new_ind).unstack('sample')
+
+del work, to_xarray, y_true, y_pred, x_true, state_ps, tmp
+
+
+
+
+
+
+

1) unit conversion / area weighting / vertical weighting#

+
+
+
# [1] Weight vertical levels by dp/g that is equivalent to a mass of air within a grid cell per unit area [kg/m2]
+# [2] Weight horizontal area of each grid cell by a[x]/mean(a[x]).
+# [3] Unit conversion to a common energy unit
+
+DS_ENERGY = {}
+for kds in ['true','pred']:
+    # Make a copy to keep original dataset
+    DS_ENERGY[kds] = DS[kds].copy(deep=True)
+
+    # vertical weighting / area weighting / unit conversion
+    for kvar in vars_mlo:
+
+        # [1] weight vertical levels by dp/g
+        #     ONLY for vertically-resolved variables, e.g., ptend_{t,q0001}
+        # dp/g = - \rho * dz
+        if vars_mlo_len[kvar] == 60:
+            DS_ENERGY[kds][kvar] = DS_ENERGY[kds][kvar] * DS_ENERGY[kds]['dp']/grav
+
+        # [2] weight area
+        #     for ALL variables
+        DS_ENERGY[kds][kvar] = DS_ENERGY[kds]['area_wgt'] * DS_ENERGY[kds][kvar]
+
+        # [3] convert units to W/m2
+        #     for variables with different units, e.g., ptend_{t,q0001}, precsc, precc
+        DS_ENERGY[kds][kvar] =  vars_mlo_energy_conv[kvar] * DS_ENERGY[kds][kvar]
+
+
+
+
+
+
+

2) Calculate metrics based on globally-averaged model output#

+
+
+
all_metrics = ['MAE','RMSE','R2']
+
+
+
+
+
+
+
# A. Calculate metrics
+# After this step,
+# ptend_{t,q0001} have [ncol, lev] dimension;
+# and the rest variables have [ncol] dimension.
+
+# if spatial analysis is desired (e.g., R2 distribution on global map or on latitude-level plane),
+# the metrics at this step should be used.
+
+
+# Select only ML output varibles
+DS_ENERGY[kds] = DS_ENERGY[kds][vars_mlo]
+
+# Caclulate 3 metrics
+Metrics = {}
+Metrics['MAE']  = (np.abs(DS_ENERGY['true']   - DS_ENERGY['pred'])).mean('time')
+Metrics['RMSE'] = np.sqrt(((DS_ENERGY['true'] - DS_ENERGY['pred'])**2.).mean('time'))
+Metrics['R2'] = 1 - ((DS_ENERGY['true'] - DS_ENERGY['pred']                    )**2.).sum('time')/\
+                    ((DS_ENERGY['true'] - DS_ENERGY['true'].mean('time'))**2.).sum('time')
+
+
+
+
+
+
+
# B. Make horizontal mean.
+# After this step,
+# ptend_{t,q0001} have [lev] dimension;
+# and the rest variables have zero dimensions, i.e., scalars.
+
+for kmetric in all_metrics:
+    Metrics[kmetric] = Metrics[kmetric].mean('ncol') # simple mean
+
+
+
+
+
+
+
# C-1. Save the result after B.
+# to save in a table format as a csv file, the level dimensions are flattened.
+
+Metrics_stacked = {}
+for kmetric in all_metrics:
+    Metrics_stacked[kmetric] = Metrics[kmetric].to_stacked_array('ml_out_idx', sample_dims='', name=kmetric)
+
+
+# save the output
+work = pd.DataFrame({'MAE':  Metrics_stacked['MAE'].values,
+                     'RMSE': Metrics_stacked['RMSE'].values,
+                     'R2':   Metrics_stacked['R2'].values}
+                    )
+work.index.name = 'output_idx'
+
+# fn_save_metrics = f'./metrics/{model_name}.metrics.csv'
+work.to_csv(fn_save_metrics)
+work
+
+
+
+
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
MAERMSER2
output_idx
00.0024880.0052180.897846
10.0049970.0099400.934334
20.0087680.0136750.966718
30.0147560.0210970.982648
40.0230480.0305240.993642
............
12334.33306572.876357-38.685671
1247.97080517.4033150.961380
12510.29917821.9544370.947612
1264.5331419.4203060.955713
1274.80630710.1187020.866159
+

128 rows × 3 columns

+
+
+
+
+
# C-2. Save the result after vertical averaging.
+# After this step,
+# ptend_{t,q0001} also have zero dimensions, i.e., scalars;
+
+# Then, the results are saved to a csv file.
+# This csv file will be used for generating plots.
+
+Metrics_vert_avg = {}
+for kmetric in all_metrics:
+    Metrics_vert_avg[kmetric] = Metrics[kmetric].mean('lev')
+    Metrics_vert_avg[kmetric] = Metrics_vert_avg[kmetric].mean('ilev') # remove dummy dim
+
+# save the output
+work = pd.DataFrame({'MAE':  Metrics_vert_avg['MAE'].to_pandas(),
+                     'RMSE': Metrics_vert_avg['RMSE'].to_pandas(),
+                     'R2':   Metrics_vert_avg['R2'].to_pandas()}
+                    )
+work.index.name = 'Variable'
+    
+# fn_save_metrics_avg = f'./metrics/{model_name}.metrics.lev-avg.csv'
+work.to_csv(fn_save_metrics_avg)
+work
+
+
+
+
+
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
MAERMSER2
Variable
ptend_t2.6827654.4209040.588987
ptend_q00014.4947527.321821-inf
cam_out_NETSW13.36091526.7075510.982684
cam_out_FLWDS5.2244686.9686420.924299
cam_out_PRECSC2.6839154.733680-inf
cam_out_PRECC34.33306572.876357-38.685671
cam_out_SOLS7.97080517.4033150.961380
cam_out_SOLL10.29917821.9544370.947612
cam_out_SOLSD4.5331419.4203060.955713
cam_out_SOLLD4.80630710.1187020.866159
+
+
+
+
+
+

Check the performance variance across top 160 MLP models (out of ~8,000 tuning trials)#

+
+
+
# vertical avg
+METRIC_MLP160 = {}
+for kmetric in ['MAE', 'RMSE', 'R2']:
+    METRIC_MLP160[kmetric] = {}
+    for kmodel in ['MLPv1', 'MLPv2']:
+        csv_path  = f'./metrics_{kmodel}-top160/*.metrics.lev-avg.csv'
+                    # this path contains 160 csv files for top 160 MLP models
+        csv_files = sorted(glob.glob(csv_path))
+        df_list = (pd.read_csv(kf, index_col=0)[kmetric] for kf in csv_files)
+        METRIC_MLP160[kmetric][kmodel] = pd.concat(df_list,axis=1)
+        
+METRIC_MLP160_STAT = {}
+for kmetric in ['MAE', 'RMSE', 'R2']:
+    METRIC_MLP160_STAT[kmetric] = {}
+    for kmodel in ['MLPv1', 'MLPv2']:
+        df_list = {'mean': METRIC_MLP160[kmetric][kmodel].mean(axis=1),
+                   'median': METRIC_MLP160[kmetric][kmodel].median(axis=1),
+                   'std': METRIC_MLP160[kmetric][kmodel].std(axis=1),
+                   'p05': METRIC_MLP160[kmetric][kmodel].quantile(0.05,axis=1),
+                   'p10': METRIC_MLP160[kmetric][kmodel].quantile(0.10,axis=1),
+                   'p25': METRIC_MLP160[kmetric][kmodel].quantile(0.25,axis=1),
+                   'p75': METRIC_MLP160[kmetric][kmodel].quantile(0.75,axis=1),
+                   'p90': METRIC_MLP160[kmetric][kmodel].quantile(0.90,axis=1),
+                   'p95': METRIC_MLP160[kmetric][kmodel].quantile(0.95,axis=1),
+                   'min': METRIC_MLP160[kmetric][kmodel].min(axis=1),
+                   'max': METRIC_MLP160[kmetric][kmodel].max(axis=1),
+                  }
+        METRIC_MLP160_STAT[kmetric][kmodel] = pd.concat(df_list,axis=1)
+        
+# Warning message will generate due to nans and infs in R2.
+
+
+
+
+
/global/homes/s/sungduk/.conda/envs/tf_conda/lib/python3.10/site-packages/numpy/lib/function_base.py:4573: RuntimeWarning: invalid value encountered in subtract
+  diff_b_a = subtract(b, a)
+/global/homes/s/sungduk/.conda/envs/tf_conda/lib/python3.10/site-packages/numpy/lib/function_base.py:4573: RuntimeWarning: invalid value encountered in subtract
+  diff_b_a = subtract(b, a)
+
+
+
+
+
+
+
# without vertical avg
+METRIC_MLP160_LEV = {}
+for kmetric in ['MAE', 'RMSE', 'R2']:
+    METRIC_MLP160_LEV[kmetric] = {}
+    for kmodel in ['MLPv1', 'MLPv2']:
+        csv_path  = f'./metrics_{kmodel}-top160/*.metrics.csv'
+                    # this path contains 160 csv files for top 160 MLP models
+        csv_files = sorted(glob.glob(csv_path))
+        df_list = (pd.read_csv(kf, index_col=0)[kmetric] for kf in csv_files)
+        METRIC_MLP160_LEV[kmetric][kmodel] = pd.concat(df_list,axis=1)
+        
+METRIC_MLP160_STAT_LEV = {}
+for kmetric in ['MAE', 'RMSE', 'R2']:
+    METRIC_MLP160_STAT_LEV[kmetric] = {}
+    for kmodel in ['MLPv1', 'MLPv2']:
+        df_list = {'mean': METRIC_MLP160_LEV[kmetric][kmodel].mean(axis=1),
+                   'median': METRIC_MLP160_LEV[kmetric][kmodel].median(axis=1),
+                   'std': METRIC_MLP160_LEV[kmetric][kmodel].std(axis=1),
+                   'p05': METRIC_MLP160_LEV[kmetric][kmodel].quantile(0.05,axis=1),
+                   'p10': METRIC_MLP160_LEV[kmetric][kmodel].quantile(0.10,axis=1),
+                   'p25': METRIC_MLP160_LEV[kmetric][kmodel].quantile(0.25,axis=1),
+                   'p75': METRIC_MLP160_LEV[kmetric][kmodel].quantile(0.75,axis=1),
+                   'p90': METRIC_MLP160_LEV[kmetric][kmodel].quantile(0.90,axis=1),
+                   'p95': METRIC_MLP160_LEV[kmetric][kmodel].quantile(0.95,axis=1),
+                   'min': METRIC_MLP160_LEV[kmetric][kmodel].min(axis=1),
+                   'max': METRIC_MLP160_LEV[kmetric][kmodel].max(axis=1),
+                  }
+        METRIC_MLP160_STAT_LEV[kmetric][kmodel] = pd.concat(df_list,axis=1)
+        
+# Warning message will generate due to nans and infs in R2.
+
+
+
+
+
/global/homes/s/sungduk/.conda/envs/tf_conda/lib/python3.10/site-packages/numpy/lib/function_base.py:4573: RuntimeWarning: invalid value encountered in subtract
+  diff_b_a = subtract(b, a)
+/global/homes/s/sungduk/.conda/envs/tf_conda/lib/python3.10/site-packages/numpy/lib/function_base.py:4573: RuntimeWarning: invalid value encountered in subtract
+  diff_b_a = subtract(b, a)
+
+
+
+
+
+
+
+

<PART 2> Load pre-saved csv files and make figures#

+
+

Plot metrics (model intercomparison)#

+

The plotting requires presaved csv files.
+One csv files per one model.

+
+
+
# file names for saved metrics (from C-2)
+fn_metrics = {'CNN':  './metrics/CNN.metrics.lev-avg.csv',
+              'MLP':  './metrics/MLP.metrics.lev-avg.csv',
+              'RPN':  './metrics/RPN.metrics.lev-avg.csv',
+              'cVAE': './metrics/cVAE.metrics.lev-avg.csv',
+              'HSR':  './metrics/HSR.metrics.lev-avg.csv',
+              'ED':  './metrics/ED.metrics.lev-avg.csv',
+             }
+
+# file names for saved metrics (from C-1)
+fn_metrics_stacked = {'CNN':  './metrics/CNN.metrics.csv',
+                      'MLP':  './metrics/MLP.metrics.csv',
+                      'RPN':  './metrics/RPN.metrics.csv',
+                      'cVAE': './metrics/cVAE.metrics.csv',
+                      'HSR':  './metrics/HSR.metrics.csv',
+                      'ED':  './metrics/ED.metrics.csv',
+                     }
+
+# assign bar colors
+# https://davidmathlogic.com/colorblind/#%23000000-%23E69F00-%2356B4E9-%23009E73-%23F0E442-%230072B2-%23D55E00-%23CC79A7
+lc_model = {'CNN':  '#0072B2',
+            'HSR':  '#E69F00',
+            'MLP':  '#2B2B2B', #'#882255',
+            'RPN':  '#009E73',
+            'cVAE': '#D55E00',
+            'ED':   '#882255' #'#56B4E9'
+           }
+
+# variable short name
+var_short_name = {'ptend_t': 'dT/dt',
+                  'ptend_q0001':'dq/dt',
+                  'cam_out_NETSW':  'NETSW',
+                  'cam_out_FLWDS':  'FLWDS',
+                  'cam_out_PRECSC': 'PRECSC',
+                  'cam_out_PRECC': 'PRECC',
+                  'cam_out_SOLS': 'SOLS',
+                  'cam_out_SOLL': 'SOLL',
+                  'cam_out_SOLSD': 'SOLSD',
+                  'cam_out_SOLLD': 'SOLLD',
+                 }
+
+# map 'output_idx' (fn_metrics_stacked) and 'variables' (fn_metrics)
+var_idx = {}
+var_idx['ptend_t']      = (0,60)
+var_idx['ptend_q0001']  = (60,120)
+var_idx['surface_vars'] = (120,128)
+
+
+
+
+
+
+
# CRPS metrics
+fn_crps = {'RPN':  './metrics_crps/RPN_CRPS.npy',
+           'cVAE': './metrics_crps/cVAE_CRPS.npy',
+           'HSR':  './metrics_crps/HSR_CRPS.npy',
+          }
+
+
+
+
+
+

plot 1. level aggregated metrics#

+
+
+

This figure is included in the main text#

+
+
+
# ordering models in a plot
+plot_this_models = ['CNN', 'ED', 'HSR','MLP','RPN','cVAE'] 
+plot_this_models_crps = ['HSR', 'RPN', 'cVAE']
+
+# which metrics to plot?
+plot_this_metrics = ['MAE', 'RMSE', 'R2', 'CRPS']
+
+# subpanel score
+abc='abcdefg'
+
+
+
+
+
+
+
# Packing metrics of different models into one Pandas dataframe
+
+PLOTDATA = {}
+for kmodel in plot_this_models:
+    # MSE, R2, RMSE
+    PLOTDATA[kmodel] = pd.read_csv(fn_metrics[kmodel], index_col=0)
+    
+    # CRPS
+    if kmodel in fn_crps.keys():
+        work = np.load(fn_crps[kmodel])
+        work = [work[slice(*var_idx['ptend_t'])].mean(),     # vertical_avg
+                work[slice(*var_idx['ptend_q0001'])].mean(), # vertical_avg
+                *work[slice(*var_idx['surface_vars'])],
+               ]
+        PLOTDATA[kmodel]['CRPS'] = work
+    else:
+        continue
+    
+
+PLOTDATA_by_METRIC = {}
+for kmetric in plot_this_metrics:
+    if kmetric in ['CRPS']:
+        k_plot_this_models = plot_this_models_crps
+    if kmetric in ['MAE', 'R2', 'RMSE']:
+        k_plot_this_models = plot_this_models
+    PLOTDATA_by_METRIC[kmetric] = pd.DataFrame([PLOTDATA[kmodel][kmetric] for kmodel in k_plot_this_models],
+                                               index=k_plot_this_models
+                                              )
+
+
+
+
+
+
+
fig, _ax = plt.subplots(nrows  = len(plot_this_metrics), 
+                        sharex = True)
+
+for k, kmetric in enumerate(plot_this_metrics):
+    ax = _ax[k]
+    plotdata = PLOTDATA_by_METRIC[kmetric]
+    plotdata = plotdata.rename(columns=var_short_name)
+    plotdata = plotdata.transpose()
+    plotdata.plot.bar(color = [lc_model[kmodel] for kmodel in plotdata.keys()],
+                      width = .4 if kmetric=='CRPS' else .6,
+                      legend = False,
+                      ax=ax)
+
+    ax.set_title(f'({abc[k]}) {kmetric}')
+    ax.set_xlabel('Output variable')
+    ax.set_xticklabels(plotdata.index, rotation=0, ha='center')
+    
+    # top-160 MLP spread
+    if True and (kmetric in METRIC_MLP160_STAT.keys()):
+        errbar_model = 'MLPv1'
+        x_off       = 0.4
+        plot_x      = np.arange(len(plotdata)) + x_off
+        plot_m      = METRIC_MLP160_STAT[kmetric][errbar_model]['median']
+        plot_lower  = plot_m - METRIC_MLP160_STAT[kmetric][errbar_model]['p05'] 
+        plot_upper  = METRIC_MLP160_STAT[kmetric][errbar_model]['p95'] - plot_m
+        
+        ax.errorbar(plot_x, plot_m, yerr=[plot_lower,plot_upper],
+                    color='k', capsize=2.5, ls='none', lw=.5)
+    
+    # no units for R2
+    if kmetric != 'R2':
+        ax.set_ylabel('W/m2')
+    
+    # not plotting negative R2 values
+    if kmetric == 'R2':
+        ax.set_ylim(0,1)
+        
+    # log y scale
+    if kmetric != 'R2':
+        ax.set_yscale('log')
+
+    fig.set_size_inches(7, 8)
+
+_ax[0].legend(ncols=3, columnspacing=.9, labelspacing=.3,
+              handleheight=.07, handlelength=1.5, handletextpad=.2,
+              borderpad=.2,
+              loc='upper right')
+
+fig.tight_layout()
+
+
+
+
+../_images/9c1f17671425f7c6b2e4305ba3542bad95fd737f7bead28feea33456a28888b5.png +
+
+
+
+

plot 2. vertical profile of metrics for tendency variables#

+
+
+

This figure is included in the SI#

+
+
+
# tell which models and which metrics you want to plot
+plot_this_models = ['CNN', 'ED', 'HSR','MLP','RPN','cVAE'] 
+plot_this_models_crps = ['HSR', 'RPN', 'cVAE']
+
+plot_this_metrics = ['MAE', 'RMSE', 'R2', 'CRPS']
+
+
+
+
+
+
+
# Packing metrics of different models into one Pandas dataframe
+
+PLOTDATA = {}
+for kmodel in plot_this_models:
+    # MSE, R2, RMSE
+    PLOTDATA[kmodel] = pd.read_csv(fn_metrics_stacked[kmodel], index_col=0)
+    
+    # CRPS
+    if kmodel in fn_crps.keys():
+        work = np.load(fn_crps[kmodel])
+        work = [*work[slice(*var_idx['ptend_t'])],
+                *work[slice(*var_idx['ptend_q0001'])],
+                *work[slice(*var_idx['surface_vars'])],
+               ]
+        PLOTDATA[kmodel]['CRPS'] = work
+    else:
+        continue
+    
+
+PLOTDATA_by_METRIC = {}
+for kmetric in plot_this_metrics:
+    if kmetric in ['CRPS']:
+        k_plot_this_models = plot_this_models_crps
+    if kmetric in ['MAE', 'R2', 'RMSE']:
+        k_plot_this_models = plot_this_models
+    PLOTDATA_by_METRIC[kmetric] = pd.DataFrame([PLOTDATA[kmodel][kmetric] for kmodel in k_plot_this_models],
+                                               index=k_plot_this_models
+                                              )
+
+
+
+
+
+
+
abc='abcdefg'
+for kvar in ['ptend_t','ptend_q0001']:
+    fig, _ax = plt.subplots(ncols=2, nrows=2)
+    _ax = _ax.flatten()
+    for k, kmetric in enumerate(plot_this_metrics):
+        ax = _ax[k]
+        idx_start = var_idx[kvar][0]
+        idx_end = var_idx[kvar][1]
+        plotdata = PLOTDATA_by_METRIC[kmetric].iloc[:,idx_start:idx_end]
+        if kvar == 'ptend_q0001':
+            plotdata.columns = plotdata.columns - 60
+        if kvar=='ptend_q0001': # this is to keep the right x axis range.
+            plotdata = plotdata.where(~np.isinf(plotdata),-999)
+        plotdata = plotdata.transpose()
+        plotdata.plot(color = [lc_model[kmodel] for kmodel in plotdata.keys()],
+                      legend=False,
+                      ax=ax,
+                     )
+        # top-160 MLP spread
+        if True and (kmetric in METRIC_MLP160_STAT_LEV.keys()):
+            errbar_model = 'MLPv1'
+            plot_x     = np.arange(0,60)
+            plot_lower = METRIC_MLP160_STAT_LEV[kmetric][errbar_model]['p05']
+            plot_lower = plot_lower.iloc[idx_start:idx_end]
+            plot_upper = METRIC_MLP160_STAT_LEV[kmetric][errbar_model]['p95']
+            plot_upper = plot_upper.iloc[idx_start:idx_end]
+            
+            ax.fill_between(plot_x, plot_lower, plot_upper,
+                            color='silver', zorder=-100)
+
+
+        ax.set_xlabel('Level index')
+        ax.set_title(f'({abc[k]}) {kmetric} ({var_short_name[kvar]})')
+        if kmetric != 'R2':
+            ax.set_ylabel('W/m2')
+
+        # R2 ylim
+        if  (kmetric=='R2'):
+            ax.set_ylim(0,1.05)
+
+    # legend
+    _ax[0].legend(ncols=1, labelspacing=.3,
+              handleheight=.07, handlelength=1.5, handletextpad=.2,
+              borderpad=.3,
+              loc='upper left')
+    
+    fig.tight_layout()
+    fig.set_size_inches(7,4.5)
+
+
+
+
+../_images/92c3a7635762776f3266559b60b1d53fe896ddf61c62009ad8f6522d32ff80a4.png +../_images/14acc70b325ebb0d917cf65415031f93f0b23e2eb18304b2cc1ec447fead3126.png +
+
+
+
+

plot 3. plot 1 and plot 2 together for MAE and R2#

+
+
+

This figure is included in the main text#

+
+
+
sw_log = False
+
+abc='abbcdefghij'
+# fig, _ax = plt.subplots(ncols=2, nrows=4,
+#                         gridspec_kw={'height_ratios': [1.6,1,1,1]})
+
+fig, _ax = plt.subplots(ncols=2, nrows=3)
+gs = _ax[0, 0].get_gridspec()
+# remove the underlying axes
+for ax in _ax[0, :]:
+    ax.remove()
+axbig = fig.add_subplot(gs[0, 0:])
+
+
+
+# top rows
+
+plot_this_models = ['CNN', 'ED', 'HSR','MLP','RPN','cVAE']
+plot_this_metrics = ['MAE']
+
+PLOTDATA = {}
+for kmodel in plot_this_models:
+    PLOTDATA[kmodel] = pd.read_csv(fn_metrics[kmodel], index_col=0)
+
+PLOTDATA_by_METRIC = {}
+for kmetric in plot_this_metrics:
+    PLOTDATA_by_METRIC[kmetric] = pd.DataFrame([PLOTDATA[kmodel][kmetric] for kmodel in plot_this_models],
+                                               index=plot_this_models
+                                              )
+
+ax = axbig
+plotdata = PLOTDATA_by_METRIC[plot_this_metrics[0]]
+plotdata = plotdata.rename(columns=var_short_name)
+plotdata = plotdata.transpose()
+plotdata.plot.bar(color=[lc_model[kmodel] for kmodel in plot_this_models],
+                  legend = False,
+                  width = 0.6,
+                  ax=ax)
+ax.set_xticklabels(plotdata.index, rotation=0, ha='center')
+ax.set_xlabel('')
+# ax.set_title(f'({abc[k]}) {kmetric}')
+ax.set_ylabel('W/m2')
+
+ax.text(0.03, 0.93, f'(a) {kmetric}', horizontalalignment='left',
+       verticalalignment='top', transform=ax.transAxes,
+       fontweight='demi')
+
+if sw_log:
+    ax.set_yscale('log')
+
+ax.legend(ncols=2,
+          columnspacing=.8,
+          labelspacing=.3,
+          handleheight=.1,
+          handlelength=1.5,
+          handletextpad=.2,
+          borderpad=.4,
+          frameon=True,
+          loc='upper right')
+
+# top-160 MLP spread
+errbar_model = 'MLPv1'
+x_off       = 0.4
+plot_x      = np.arange(len(plotdata)) + x_off
+plot_m      = METRIC_MLP160_STAT[plot_this_metrics[0]][errbar_model]['median']
+plot_lower  = plot_m - METRIC_MLP160_STAT[plot_this_metrics[0]][errbar_model]['p05'] 
+plot_upper  = METRIC_MLP160_STAT[plot_this_metrics[0]][errbar_model]['p95'] - plot_m
+
+
+ax.errorbar(plot_x, plot_m, yerr=[plot_lower,plot_upper],
+            color='k', capsize=2.5, ls='none', lw=.5)
+
+# bottom rows
+
+plot_this_models = ['CNN', 'ED', 'HSR','MLP','RPN','cVAE']
+lot_this_models = ['ED', 'HSR','MLP']
+plot_this_metrics = ['MAE', 'R2']
+
+var_idx = {}
+var_idx['ptend_t'] = (0,60)
+var_idx['ptend_q0001'] = (60,120)
+
+PLOTDATA = {}
+for kmodel in plot_this_models:
+    PLOTDATA[kmodel] = pd.read_csv(fn_metrics_stacked[kmodel], index_col=0)
+
+PLOTDATA_by_METRIC = {}
+for kmetric in plot_this_metrics:
+    PLOTDATA_by_METRIC[kmetric] = pd.DataFrame([PLOTDATA[kmodel][kmetric] for kmodel in plot_this_models],
+                                               index=plot_this_models
+                                              )
+
+
+for kk, kvar in enumerate(['ptend_t','ptend_q0001']):
+    for k, kmetric in enumerate(plot_this_metrics):
+        ax = _ax[k+1, 0 if kvar=='ptend_t' else 1]
+        idx_start = var_idx[kvar][0]
+        idx_end = var_idx[kvar][1]
+        plotdata = PLOTDATA_by_METRIC[kmetric].iloc[:,idx_start:idx_end]
+        if kvar == 'ptend_q0001':
+            plotdata.columns = plotdata.columns - 60
+        if kvar=='ptend_q0001': # this is to keep the right x axis range.
+            plotdata = plotdata.where(~np.isinf(plotdata),-999)
+        plotdata.transpose()\
+        .plot(color=[lc_model[kmodel] for kmodel in plot_this_models],
+              legend=False,
+              ax=ax,
+             )
+
+        # top-160 MLP spread
+        if True and (kmetric in METRIC_MLP160_STAT_LEV.keys()):
+            errbar_model = 'MLPv1'
+            plot_x     = np.arange(0,60)
+            plot_lower = METRIC_MLP160_STAT_LEV[kmetric][errbar_model]['p05']
+            plot_lower = plot_lower.iloc[idx_start:idx_end]
+            plot_upper = METRIC_MLP160_STAT_LEV[kmetric][errbar_model]['p95']
+            plot_upper = plot_upper.iloc[idx_start:idx_end]
+            
+            ax.fill_between(plot_x, plot_lower, plot_upper,
+                            color='silver', zorder=-100)
+
+        #ax.set_title(f'({abc[k]}) {kmetric}')
+        if k==0:
+            ax.set_ylabel(f'W/m2')
+            ax.set_xlabel('')
+            ax.set_xticklabels('')
+        elif k==1:
+            ax.set_xlabel('Level index')
+
+        
+        if abc[kk+2*k+2] == 'd':
+            ax.text(0.97, 0.93, f'({abc[kk+2*k+2]}) {kmetric}, {var_short_name[kvar]}', horizontalalignment='right',
+                   verticalalignment='top', transform=ax.transAxes,
+                   fontweight='demi')
+        else:
+            ax.text(0.03, 0.93, f'({abc[kk+2*k+2]}) {kmetric}, {var_short_name[kvar]}', horizontalalignment='left',
+                   verticalalignment='top', transform=ax.transAxes,
+                   fontweight='demi')
+
+        if sw_log:
+            ax.set_yscale('log')
+
+        # R2 ylim
+        if  (kmetric=='R2'):
+            ax.set_ylim(0,1.05)
+
+fig.set_size_inches(9,5.25)
+
+
+
+
+../_images/6ae8a68f242679e29d61ef6240d1480ecedb1319f0020a53724ea85bc79539b2.png +
+
+
+
+

plot 4. scatter plot#

+

(Make sure to run the cells upto Cell#8)

+
+
+
from matplotlib import cm, colors
+from mpl_toolkits.axes_grid1 import ImageGrid
+from sklearn.metrics import r2_score
+
+
+
+
+
+
+
# load model prediction and pack into a single dictionary.
+DS_HIST = {}
+for kmodel in ['true', 'CNN', 'ED', 'HSR', 'MLP', 'RPN', 'cVAE']:
+           # make sure that 'true' should be the first entry
+    if kmodel=='true':
+        fn_y = '../npy_data_conversion/npy_files/val_target_stride6.npy'
+        work = np.load(fn_y).astype(np.float64)
+        fn_x = '../npy_data_conversion/npy_files/val_input_stride6.npy'
+        x_true = np.load(fn_x).astype(np.float64)
+    else:
+        fn_y = FN_MODEL_OUTPUT[kmodel]
+        if fn_y[-3:] == '.h5':
+            work = xr.open_dataset(fn_y)['pred'].values
+        else:
+            work = np.load(fn_y).astype(np.float64)
+    N_samples = work.shape[0] 
+
+    # [1] Construct dictionary for xarray dataset
+    #     format is key for variable name /
+    #               value for a turple of (dimension names, data).
+    to_xarray = {}
+    for k, kvar in enumerate(vars_mlo):
+
+        # length of variable (ie, number of levels)
+        kvar_len = vars_mlo_len[kvar]
+
+        # set dimensions of variable
+        if kvar_len == 60:
+            kvar_dims = (dim_name_sample, dim_name_level)
+        elif kvar_len == 1:
+            kvar_dims = dim_name_sample
+
+        # set start and end indices of variable in the loaded numpy array
+        # then, add 'kvar':(kvar_dims, <np_array>) to dictionary
+        if k==0: ind1=0
+        ind2 = ind1 + kvar_len
+
+        # scaled output
+        kvar_data = np.squeeze(work[:,ind1:ind2])
+        # # unscaled output
+        # kvar_data = kvar_data / mlo_scale[kvar].values
+
+        to_xarray[kvar] = (kvar_dims, kvar_data)
+
+        ind1 = ind2
+
+    # [2] convert dict to xarray dataset
+    DS_HIST[kmodel] = xr.Dataset(to_xarray)
+
+    # [3] add surface pressure ('state_ps') from ml input
+    # normalized ps
+    state_ps =  xr.DataArray(x_true[:,120], dims=('sample'), name='state_ps')
+    # denormalized ps
+    state_ps = state_ps * (mli_max['state_ps'] - mli_min['state_ps']) + mli_mean['state_ps']
+    DS_HIST[kmodel]['state_ps'] = state_ps
+
+    # [4] add grid information
+    DS_HIST[kmodel] = xr.merge([DS_HIST[kmodel], ds_grid])
+
+    # [5] add pressure thickness of each level, dp
+    # FYI, in a hybrid sigma vertical coordinate system, pressure at level z is
+    # P[x,z] = hyam[z]*P0 + hybm[z]*PS[x,z],
+    # where, hyam and hybm are 
+    tmp = DS_HIST[kmodel]['P0']*DS_HIST[kmodel]['hyai'] \
+          + DS_HIST[kmodel]['state_ps']*DS_HIST[kmodel]['hybi']
+    tmp = tmp.isel(ilev=slice(1,61)).values - tmp.isel(ilev=slice(0,60)).values
+    tmp = tmp.transpose()
+    DS_HIST[kmodel]['dp'] = xr.DataArray(tmp, dims=('sample', 'lev'))
+
+    # [6] break (sample) to (ncol,time)
+    N_timestep = int(N_samples/N_ncol)
+    dim_ncol     = np.arange(N_ncol)
+    dim_timestep = np.arange(N_timestep)
+    new_ind = pd.MultiIndex.from_product([dim_timestep, dim_ncol],
+                                         names=['time', 'ncol'])
+    DS_HIST[kmodel] = DS_HIST[kmodel].assign_coords(sample=new_ind).unstack('sample')
+
+del work, to_xarray, x_true, state_ps, tmp
+
+
+
+
+
+
+
## PLOT ##
+def plot_hexbin(kvar:str, klev=-1, 
+                xmin=0., xmax=3., ymin=0., ymax=3.):
+    fig = plt.figure()
+    imgrid = ImageGrid(fig, 111,  # similar to subplot(111)
+                       nrows_ncols=(1, 6),  # creates 2x2 grid of axes
+                       axes_pad=0.2,  # pad between axes in inch.
+                       cbar_mode='single',
+                       cbar_location='right',
+                       cbar_pad=0.2,
+                       aspect=True
+                       )
+
+    panel_id='abcdefg'
+    plot_this_models = ['CNN', 'ED', 'HSR','MLP','RPN','cVAE'] 
+    for k, kmodel in enumerate(plot_this_models):
+
+        # read data
+        plotx = DS_HIST['true'][kvar]
+        ploty = DS_HIST[kmodel][kvar]
+        
+        if klev >= 0:
+            plotx = plotx.isel(lev=klev)
+            ploty = ploty.isel(lev=klev)
+            
+        plotx = plotx.values
+        ploty = ploty.values
+        
+        # metric calculation
+        plotrmse = np.sqrt(((plotx - ploty)**2.).mean())
+        plotR2 = r2_score(plotx, ploty, sample_weight=None, multioutput='uniform_average')
+
+        # # normalize
+        # plotx = plotx / plotx.max()
+        # ploty = ploty / plotx.max()
+
+        # n_samples
+        n      = DS_HIST['true'][kvar].values.size
+
+        ax = imgrid[k]
+        h = ax.hexbin(plotx, ploty, 
+                      gridsize=50, extent=(xmin,xmax,ymin,ymax),
+                      bins='log', 
+                      mincnt=1, vmin=1, vmax=n/125, 
+                      cmap='viridis', alpha=1.)
+        ax.axline([0,0], slope=1, color='k', linewidth=1)
+        ax.set_facecolor('#EFEFEF')
+
+
+        ax.set_xlim(xmin, xmax)
+        ax.set_ylim(ymin, ymax)
+        
+        ax.set_ylabel('Predicted')
+        ax.set_xlabel('True')
+
+        props = dict(boxstyle='round', facecolor='wheat', alpha=.7)
+        # ax.text(.05,.95, f'({panel_id[k]}) {kmodel}\nR2:{PLOTDATA_by_METRIC["R2"].loc[kmodel,kvar]:.2f}\nRMSE:{PLOTDATA_by_METRIC["RMSE"].loc[kmodel,kvar]:.2f}', 
+        #         ha='left', va='top', transform=ax.transAxes, bbox=props, fontsize=9)
+        ax.text(.05,.95, f'({panel_id[k]}) {kmodel}\nRMSE:{plotrmse:4.2E}\nR2:{plotR2:.3f}', 
+                ha='left', va='top', transform=ax.transAxes, bbox=props, fontsize=6.5)
+
+    cb = imgrid.cbar_axes[0].colorbar(h)
+    cb.set_label('Counts')
+    
+    fig.suptitle(vars_longname[kvar] + (f' (level={klev})' if klev>=0 else ''),
+                 y=.95)
+
+    #fig.set_size_inches(9,4.55) # nrows_ncols=(2, 3)
+    fig.set_size_inches(10,1.8) # nrows_ncols=(2, 3)
+    fig.set_facecolor('w')
+        
+    return
+
+
+
+
+
+
+
# plot 2d variables
+for kvar in ['cam_out_NETSW', 'cam_out_FLWDS', 'cam_out_PRECSC', 'cam_out_PRECC', 'cam_out_SOLS', 'cam_out_SOLL', 'cam_out_SOLSD', 'cam_out_SOLLD']:
+    plot_hexbin(kvar,
+            xmin=0, xmax=2.8,
+            ymin=0, ymax=2.8)
+
+
+
+
+../_images/5ea863f8636c3e4dd01bb801f3e516b85bf710b348131f4fed3564b26eaa678c.png +../_images/2ee4a8b25aa9c415fcb6ca4b3ce87fc4fff6a628e8590f64e0eb625c564ba8d5.png +../_images/01e9ea55517d61d3e122adf8066f63ec0c7de87f6b892997dfee23da30164393.png +../_images/8ee0384f380a67ce378746757e4de5e6ee69d1e05bd12634246d950cf35ac82a.png +../_images/694ba9a105e04acaaf2b0e4def9117cc6fd8f138e632c99299bc5057794b96de.png +../_images/af3470be5bc7dfd72baff121ef9bb979b14b8d6793636d374f2eec0e3950dc74.png +../_images/a591cae313c30e179519233942f128da6d95382cfabb1036e6555ae3641c75bc.png +../_images/53159bd2d16373b3a3fa7cf3d6bdc12bccb5439744acf64278a7102eeadbcff2.png +
+
+
+
+
# plot 3d variables
+plot_hexbin('ptend_t', klev=5,
+            xmin=-.35, xmax=.35,
+            ymin=-.35, ymax=.35)
+plot_hexbin('ptend_t', klev=15,
+            xmin=-.35, xmax=.35,
+            ymin=-.35, ymax=.35)
+plot_hexbin('ptend_t', klev=35,
+            xmin=-1.5, xmax=1.5,
+            ymin=-1.5, ymax=1.5)
+plot_hexbin('ptend_t', klev=55,
+            xmin=-1.5, xmax=1.5,
+            ymin=-1.5, ymax=1.5)
+
+plot_hexbin('ptend_q0001', klev=5,
+            xmin=-.0025, xmax=.0025,
+            ymin=-.0025, ymax=.0025)
+plot_hexbin('ptend_q0001', klev=15,
+            xmin=-.0065, xmax=.0065,
+            ymin=-.0065, ymax=.0065)
+plot_hexbin('ptend_q0001', klev=35,
+            xmin=-1.5, xmax=1.5,
+            ymin=-1.5, ymax=1.5)
+plot_hexbin('ptend_q0001', klev=55,
+            xmin=-2.75, xmax=2.75,
+            ymin=-2.75, ymax=2.75)
+
+
+
+
+../_images/3cd02363e952465d5d10ced24b44438ccf99753d91f7fa99d42dc7d2b0922417.png +../_images/d48027406e135b4796a5eb35b752b490270d090b6979700c442d50862585a41b.png +../_images/7e20b4b3eaa6402e3bfa52aed9713bf76e58496d12ba64cbf23712a5a83b2458.png +../_images/573f9adc1b34ea7834aa73140a19734575c688066d2948474d0653e45277a570.png +../_images/4a931fc17178b90d47872e223f82eba9325ad9b76f14f8a24ea10f8e64148ce2.png +../_images/d5afbd437e3020dc731c87756606ddc4d1ce4537eb8a9f56f3f1fda95fb25c5b.png +../_images/b8e6913d88374b49781ecdc2443fa9fb737d0a8689b83818e779eea0cc2883ca.png +../_images/1e916c961bc0efbcff91c308fb52058960bcd0256dbebdc82f3f2ea0ceebf06d.png +
+
+
+
+
+ + + + +
+ + + + +
+ +
+ +
+
+
+ +
+ +
+ +
+ + + +
+ +
+
+ Contents +
+
+ +
+ + +
+
+ +
+ +
+ +

+By The ClimSim Maintainers +

+ +
+ +
+ +

+ + © Copyright 2023. +
+ +

+ +
+ +
+ +
+ +
+ +
+ +
+
+ + +
+
+
+ + + + + + + + \ No newline at end of file diff --git a/evaluation/bias_calculation.html b/evaluation/bias_calculation.html new file mode 100644 index 0000000..2d58677 --- /dev/null +++ b/evaluation/bias_calculation.html @@ -0,0 +1,670 @@ + + + + + + + + + + + + Import data_utils — ClimSim + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+
+
+ + + Ctrl+K +
+
+ + + +
+
+ +
+ + + +
+ + + + +
+ +
+
+ + +
+
+
+ + +
+
+ + + + +
+ +
+ + + +
+ +
+
+ +
+
+ +
+ +
+ +
+ + +
+ +
+ +
+ + + + + +
+ + +
+ + + + + + + + + + + +
+ +
+ +
+
+ + + +
+

Import data_utils

+ +
+ +
+
+ + + + +
+ +
+

Import data_utils#

+
+
+
from climsim_utils.data_utils import *
+
+
+
+
+
+
+

Instantiate class#

+
+
+
grid_path = '/ocean/projects/atm200007p/jlin96/neurips_proj/ClimSim/grid_info/ClimSim_low-res_grid-info.nc'
+norm_path = '/ocean/projects/atm200007p/jlin96/neurips_proj/ClimSim/preprocessing/normalizations/'
+
+grid_info = xr.open_dataset(grid_path)
+input_mean = xr.open_dataset(norm_path + 'inputs/input_mean.nc')
+input_max = xr.open_dataset(norm_path + 'inputs/input_max.nc')
+input_min = xr.open_dataset(norm_path + 'inputs/input_min.nc')
+output_scale = xr.open_dataset(norm_path + 'outputs/output_scale.nc')
+
+data = data_utils(grid_info = grid_info, 
+                  input_mean = input_mean, 
+                  input_max = input_max, 
+                  input_min = input_min, 
+                  output_scale = output_scale)
+
+
+
+
+
+
+

Load data and set pressure grid#

+
+
+
# paths to scoring data
+input_path = '/ocean/projects/atm200007p/jlin96/neurips_proj/e3sm_train_npy/scoring_input.npy'
+target_path = '/ocean/projects/atm200007p/jlin96/neurips_proj/e3sm_train_npy/scoring_target.npy'
+
+# paths to model predictions
+cvae_pred_path = '/ocean/projects/atm200007p/shared/neurips_proj/final_metrics/predictions/cVAE/cvae_preds_manual.h5'
+ed_pred_path = '/ocean/projects/atm200007p/behrens/ED_Behrens_2022/ED_ClimSIM_1_3_pred.npy'
+hsr_pred_path = '/ocean/projects/atm200007p/shared/neurips_proj/final_metrics/predictions/HSR/hsr_preds_bestcrps.h5'
+rpn_pred_path = '/ocean/projects/atm200007p/jlin96/neurips_proj/figure_ingredients/rpn_pred_v1_stride6.npy'
+cnn_pred_path = '/ocean/projects/atm200007p/jlin96/neurips_proj/figure_ingredients/val_predict_cnn_reshaped_stride6_FINAL.npy'
+mlp_pred_path = '/ocean/projects/atm200007p/jlin96/neurips_proj/figure_ingredients/001_backup_phase-7_retrained_models_step2_lot-147_trial_0027.best.h5.npy'
+
+# set variables to V1 subset
+data.set_to_v1_vars()
+
+# path to target input
+data.input_scoring = np.load(input_path)
+
+# path to target output
+data.target_scoring = np.load(target_path)
+
+# set pressure weights
+data.set_pressure_grid(input_arr = data.input_scoring)
+
+# load model predictions
+data.model_names = ['CNN','cVAE','ED','HSR','MLP', 'RPN']
+preds = [data.load_npy_file(load_path = cnn_pred_path), 
+         data.load_h5_file(load_path = cvae_pred_path),
+         data.load_npy_file(load_path = ed_pred_path), 
+         data.load_h5_file(load_path = hsr_pred_path), 
+         data.load_npy_file(load_path = mlp_pred_path), 
+         data.load_npy_file(load_path = rpn_pred_path)]
+data.preds_scoring = dict(zip(data.model_names, preds))
+
+
+
+
+
+
+

Weight outputs#

+
    +

  1. Undo output scaling

    2. +

  2. Weight vertical levels by dp/g

    3. +

  3. Weight horizontal area of each grid cell by a[x]/mean(a[x])

    4. +

  4. Convert units to a common energy unit

    5. +
+
+
+
data.reweight_target(data_split = 'scoring')
+data.reweight_preds(data_split = 'scoring')
+
+
+
+
+
+
+

Calculate biases#

+
+
+
results_complete = {}
+for model in data.model_names:
+    results_temp = {}
+    for var in data.target_vars:
+        results_temp[var] = data.calc_bias(data.preds_weighted_scoring[model][var], 
+                                           data.target_weighted_scoring[var],
+                                           avg_grid = False)
+    results_complete[model] = results_temp
+
+
+
+
+
+
+

Get pressure axis#

+
+
+
# get plotting pressure grid
+data.data_path = '/ocean/projects/atm200007p/jlin96/neurips_proj/e3sm_train/'
+
+# set regular expressions for selecting data
+data.set_regexps(data_split = 'scoring', regexps = ['E3SM-MMF.mli.0008-0[23456789]-*-*.nc', 
+                                                    'E3SM-MMF.mli.0008-1[012]-*-*.nc', 
+                                                    'E3SM-MMF.mli.0009-01-*-*.nc'])
+# set temporal subsampling
+data.set_stride_sample(data_split = 'scoring', stride_sample = 6)
+
+# create list of files to extract data from
+data.set_filelist(data_split = 'scoring')
+
+# create pressure grid
+pressure_grid_plotting = data.get_pressure_grid_plotting(data_split = 'scoring')
+
+# average pressure grid across grid
+y_pressure = pressure_grid_plotting.mean(axis = 1)/100
+
+
+
+
+
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 4380/4380 [00:32<00:00, 135.51it/s]
+
+
+
+
+
+
+

Make plots#

+
+
+
%config InlineBackend.figure_format = 'retina'
+letters = string.ascii_lowercase
+
+
+
+
+
+
+
plt.plot(results_complete['CNN']['ptend_t'].mean(axis = 0), y_pressure)
+plt.plot(results_complete['ED']['ptend_t'].mean(axis = 0), y_pressure)
+plt.plot(results_complete['HSR']['ptend_t'].mean(axis = 0), y_pressure)
+plt.plot(results_complete['MLP']['ptend_t'].mean(axis = 0), y_pressure)
+plt.plot(results_complete['RPN']['ptend_t'].mean(axis = 0), y_pressure)
+plt.plot(results_complete['cVAE']['ptend_t'].mean(axis = 0), y_pressure)
+
+plt.legend(['CNN', 'ED', 'HSR', 'MLP', 'RPN', 'cVAE'])
+plt.xlim(-1.4, 1.4)
+plt.gca().invert_yaxis()
+plt.xlabel('$W/m^2$')
+plt.ylabel('pressure (hPa)')
+plt.title('heating bias')
+
+
+
+
+
Text(0.5, 1.0, 'heating bias')
+
+
+../_images/29b498c276c50f41200b73dea8cea619bb058a6d23311f917fda21d73f6f3675.png +
+
+
+
+
plt.plot(results_complete['CNN']['ptend_q0001'].mean(axis = 0), y_pressure)
+plt.plot(results_complete['ED']['ptend_q0001'].mean(axis = 0), y_pressure)
+plt.plot(results_complete['HSR']['ptend_q0001'].mean(axis = 0), y_pressure)
+plt.plot(results_complete['MLP']['ptend_q0001'].mean(axis = 0), y_pressure)
+plt.plot(results_complete['RPN']['ptend_q0001'].mean(axis = 0), y_pressure)
+plt.plot(results_complete['cVAE']['ptend_q0001'].mean(axis = 0), y_pressure)
+
+plt.legend(['CNN', 'ED', 'HSR', 'MLP', 'RPN', 'cVAE'])
+plt.xlim(-1.4, 1.4)
+plt.gca().invert_yaxis()
+plt.xlabel('$W/m^2$')
+plt.ylabel('pressure (hPa)')
+plt.title('moistening bias')
+
+
+
+
+
Text(0.5, 1.0, 'moistening bias')
+
+
+../_images/123ba4b98a7b14473ad36bf0407ba63b2fd1b18268dbec942a0ed52fd702528e.png +
+
+
+
+
data.metrics_names = ['bias']
+data.create_metrics_df(data_split = 'scoring')
+bias_df = pd.DataFrame([data.metrics_var_scoring[model]['bias'] for model in data.model_names],
+                        index=data.model_names)
+bias_df = bias_df.rename(columns = data.var_short_names).transpose()
+
+
+
+
+
+
+
bias_df.plot.bar()
+plt.xlabel('Output variable')
+plt.ylabel('$W/m^2$')
+plt.title('Bias')
+
+
+
+
+
Text(0.5, 1.0, 'Bias')
+
+
+../_images/68073e0f59727d197ee21594196d09a7b932120eab6a057caf030ba455fdd12f.png +
+
+
+ + + + +
+ + + + +
+ +
+ +
+
+
+ +
+ +
+ +
+ + + +
+ +
+
+ Contents +
+
+ +
+ + +
+
+ +
+ +
+ +

+By The ClimSim Maintainers +

+ +
+ +
+ +

+ + © Copyright 2023. +
+ +

+ +
+ +
+ +
+ +
+ +
+ +
+
+ + +
+
+
+ + + + + + + + \ No newline at end of file diff --git a/evaluation/main_figure_generation.html b/evaluation/main_figure_generation.html new file mode 100644 index 0000000..d2aecea --- /dev/null +++ b/evaluation/main_figure_generation.html @@ -0,0 +1,858 @@ + + + + + + + + + + + + Evaluation metrics — ClimSim + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+
+
+ + + Ctrl+K +
+
+ + + +
+
+ +
+ + + +
+ + + + +
+ +
+
+ + +
+
+
+ + +
+
+ + + + +
+ +
+ + + +
+ +
+
+ +
+
+ +
+ +
+ +
+ + +
+ +
+ +
+ + + + + +
+ + +
+ + + + + + + + + + + +
+ +
+ +
+
+ + + +
+

Evaluation metrics

+ +
+ +
+
+ + + + +
+ +
+

Evaluation metrics#

+
+

Import data_utils#

+
+
+
from climsim_utils.data_utils import *
+
+
+
+
+
2023-08-21 14:36:29.941637: I tensorflow/core/platform/cpu_feature_guard.cc:182] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations.
+To enable the following instructions: AVX2 FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags.
+2023-08-21 14:36:30.922437: W tensorflow/compiler/tf2tensorrt/utils/py_utils.cc:38] TF-TRT Warning: Could not find TensorRT
+
+
+
+
+
+
+

Instantiate class#

+
+
+
grid_path = '/ocean/projects/atm200007p/jlin96/neurips_proj/ClimSim/grid_info/ClimSim_low-res_grid-info.nc'
+norm_path = '/ocean/projects/atm200007p/jlin96/neurips_proj/ClimSim/normalizations/'
+
+grid_info = xr.open_dataset(grid_path)
+input_mean = xr.open_dataset(norm_path + 'inputs/input_mean.nc')
+input_max = xr.open_dataset(norm_path + 'inputs/input_max.nc')
+input_min = xr.open_dataset(norm_path + 'inputs/input_min.nc')
+output_scale = xr.open_dataset(norm_path + 'outputs/output_scale.nc')
+
+data = data_utils(grid_info = grid_info, 
+                  input_mean = input_mean, 
+                  input_max = input_max, 
+                  input_min = input_min, 
+                  output_scale = output_scale)
+
+
+
+
+
+
+

Load data and set pressure grid#

+
+
+
# paths to scoring data
+input_path = '/ocean/projects/atm200007p/jlin96/neurips_proj/e3sm_train_npy/scoring_input.npy'
+target_path = '/ocean/projects/atm200007p/jlin96/neurips_proj/e3sm_train_npy/scoring_target.npy'
+
+# paths to model predictions
+cvae_pred_path = '/ocean/projects/atm200007p/shared/neurips_proj/final_metrics/predictions/cVAE/cvae_preds_manual.h5'
+ed_pred_path = '/ocean/projects/atm200007p/behrens/ED_Behrens_2022/ED_ClimSIM_1_3_pred.npy'
+hsr_pred_path = '/ocean/projects/atm200007p/shared/neurips_proj/final_metrics/predictions/HSR/hsr_preds_bestcrps.h5'
+rpn_pred_path = '/ocean/projects/atm200007p/jlin96/neurips_proj/figure_ingredients/rpn_pred_v1_stride6.npy'
+cnn_pred_path = '/ocean/projects/atm200007p/jlin96/neurips_proj/figure_ingredients/val_predict_cnn_reshaped_stride6_FINAL.npy'
+mlp_pred_path = '/ocean/projects/atm200007p/jlin96/neurips_proj/figure_ingredients/001_backup_phase-7_retrained_models_step2_lot-147_trial_0027.best.h5.npy'
+
+# set variables to V1 subset
+data.set_to_v1_vars()
+
+# path to target input
+data.input_scoring = np.load(input_path)
+
+# path to target output
+data.target_scoring = np.load(target_path)
+
+# set pressure weights
+data.set_pressure_grid(input_arr = data.input_scoring)
+
+# load model predictions
+data.model_names = ['CNN','cVAE','ED','HSR','MLP', 'RPN']
+preds = [data.load_npy_file(load_path = cnn_pred_path), 
+         data.load_h5_file(load_path = cvae_pred_path),
+         data.load_npy_file(load_path = ed_pred_path), 
+         data.load_h5_file(load_path = hsr_pred_path), 
+         data.load_npy_file(load_path = mlp_pred_path), 
+         data.load_npy_file(load_path = rpn_pred_path)]
+data.preds_scoring = dict(zip(data.model_names, preds))
+
+
+
+
+
+
+

Weight outputs#

+
    +

  1. Undo output scaling

    2. +

  2. Weight vertical levels by dp/g

    3. +

  3. Weight horizontal area of each grid cell by a[x]/mean(a[x])

    4. +

  4. Convert units to a common energy unit

    5. +
+
+
+
data.reweight_target(data_split = 'scoring')
+data.reweight_preds(data_split = 'scoring')
+
+data.metrics_names = ['MAE', 'RMSE', 'R2']
+data.create_metrics_df(data_split = 'scoring')
+
+
+
+
+
/ocean/projects/atm200007p/jlin96/neurips_proj/ClimSim/climsim_utils/data_utils.py:649: RuntimeWarning: divide by zero encountered in divide
+  r_squared = 1 - sq_diff.sum(axis = 0)/tss_time.sum(axis = 0) # sum over time
+/ocean/projects/atm200007p/jlin96/neurips_proj/ClimSim/climsim_utils/data_utils.py:649: RuntimeWarning: invalid value encountered in divide
+  r_squared = 1 - sq_diff.sum(axis = 0)/tss_time.sum(axis = 0) # sum over time
+/ocean/projects/atm200007p/jlin96/neurips_proj/ClimSim/climsim_utils/data_utils.py:649: RuntimeWarning: divide by zero encountered in divide
+  r_squared = 1 - sq_diff.sum(axis = 0)/tss_time.sum(axis = 0) # sum over time
+/ocean/projects/atm200007p/jlin96/neurips_proj/ClimSim/climsim_utils/data_utils.py:649: RuntimeWarning: divide by zero encountered in divide
+  r_squared = 1 - sq_diff.sum(axis = 0)/tss_time.sum(axis = 0) # sum over time
+/ocean/projects/atm200007p/jlin96/neurips_proj/ClimSim/climsim_utils/data_utils.py:649: RuntimeWarning: divide by zero encountered in divide
+  r_squared = 1 - sq_diff.sum(axis = 0)/tss_time.sum(axis = 0) # sum over time
+/ocean/projects/atm200007p/jlin96/neurips_proj/ClimSim/climsim_utils/data_utils.py:649: RuntimeWarning: divide by zero encountered in divide
+  r_squared = 1 - sq_diff.sum(axis = 0)/tss_time.sum(axis = 0) # sum over time
+/ocean/projects/atm200007p/jlin96/neurips_proj/ClimSim/climsim_utils/data_utils.py:649: RuntimeWarning: invalid value encountered in divide
+  r_squared = 1 - sq_diff.sum(axis = 0)/tss_time.sum(axis = 0) # sum over time
+/ocean/projects/atm200007p/jlin96/neurips_proj/ClimSim/climsim_utils/data_utils.py:649: RuntimeWarning: divide by zero encountered in divide
+  r_squared = 1 - sq_diff.sum(axis = 0)/tss_time.sum(axis = 0) # sum over time
+
+
+
+
+
+
+

Create plots#

+
+
+
%config InlineBackend.figure_format = 'retina'
+letters = string.ascii_lowercase
+
+
+
+
+
+
+
dict_idx = data.metrics_var_scoring
+dict_var = data.metrics_idx_scoring
+
+plot_df_byvar = {}
+for metric in data.metrics_names:
+    plot_df_byvar[metric] = pd.DataFrame([dict_idx[model][metric] for model in data.model_names],
+                                               index=data.model_names)
+    plot_df_byvar[metric] = plot_df_byvar[metric].rename(columns = data.var_short_names).transpose()
+
+
+
+
+
+
+
plot_df_byvar
+
+
+
+
+
{'MAE':                 CNN       cVAE         ED        HSR        MLP        RPN
+ variable                                                                  
+ $dT/dt$    2.585418   2.732051   2.864516   2.844789   2.682765   2.685253
+ $d1/dt$    4.400639   4.679666   4.673036   4.784426   4.494752   4.592410
+ NETSW     18.845029  19.729191  14.967681  19.816967  13.360914  18.878685
+ FLWDS      8.597703   6.588051   6.893821   6.266635   5.224468   6.017693
+ PRECSC     3.364180   3.321645   3.045900   3.511187   2.683915   3.328444
+ PRECC     37.825799  38.805963  37.250200  42.379195  34.333065  37.455143
+ SOLS      10.830984  10.939536   8.553734  11.314418   7.970805  10.355971
+ SOLL      13.152971  13.461330  10.923703  13.602788  10.299178  12.955313
+ SOLSD      5.816831   6.158717   5.074686   6.330680   4.533141   5.845594
+ SOLLD      5.678992   6.065650   5.136236   6.214983   4.806307   5.701871,
+ 'RMSE':                 CNN       cVAE         ED        HSR        MLP        RPN
+ variable                                                                  
+ $dT/dt$    4.369106   4.720812   4.696491   4.824681   4.420904   4.481825
+ $d1/dt$    7.283705   7.779619   7.642989   7.895550   7.321821   7.518002
+ NETSW     36.906747  38.364034  28.537017  37.773019  26.707550  33.596484
+ FLWDS     10.854747   8.529763   9.070148   8.219851   6.968642   7.913915
+ PRECSC     6.000881   6.182026   5.078154   6.094945   4.733680   5.511482
+ PRECC     85.311686  88.711142  76.681968  90.640413  72.876357  76.578238
+ SOLS      22.916411  23.264985  17.998650  23.614026  17.403315  20.606768
+ SOLL      27.251214  27.812571  22.540191  27.778323  21.954437  25.218968
+ SOLSD     12.126842  12.639114   9.916940  12.395840   9.420306  10.996472
+ SOLLD     12.102328  12.624935  10.416563  12.474274  10.118702  11.247108,
+ 'R2':                CNN      cVAE         ED        HSR        MLP        RPN
+ variable                                                                
+ $dT/dt$   0.627260  0.589693   0.541714   0.567524   0.588987   0.617206
+ $d1/dt$       -inf      -inf       -inf       -inf       -inf       -inf
+ NETSW     0.943800  0.956531   0.980463   0.959136   0.982684   0.968036
+ FLWDS     0.827836  0.883237   0.801872   0.904115   0.924299   0.912173
+ PRECSC         NaN      -inf       -inf       -inf        NaN       -inf
+ PRECC     0.076930 -0.926289 -17.908577 -68.347502 -38.685671 -67.943146
+ SOLS      0.927029  0.929375   0.959621   0.928777   0.961380   0.942723
+ SOLL      0.915685  0.915267   0.945472   0.916143   0.947612   0.928029
+ SOLSD     0.927229  0.921261   0.950717   0.922621   0.955713   0.939980
+ SOLLD     0.812537  0.796327   0.857180   0.797273   0.866159   0.836630}
+
+
+
+
+
+
+
dict_idx
+
+
+
+
+
{'CNN':                       MAE       RMSE        R2
+ variable                                      
+ ptend_t          2.585418   4.369106   0.62726
+ ptend_q0001      4.400639   7.283705      -inf
+ cam_out_NETSW   18.845029  36.906747    0.9438
+ cam_out_FLWDS    8.597703  10.854747  0.827836
+ cam_out_PRECSC    3.36418   6.000881       NaN
+ cam_out_PRECC   37.825799  85.311686   0.07693
+ cam_out_SOLS    10.830984  22.916411  0.927029
+ cam_out_SOLL    13.152971  27.251214  0.915685
+ cam_out_SOLSD    5.816831  12.126842  0.927229
+ cam_out_SOLLD    5.678992  12.102328  0.812537,
+ 'cVAE':                       MAE       RMSE        R2
+ variable                                      
+ ptend_t          2.732051   4.720812  0.589693
+ ptend_q0001      4.679666   7.779619      -inf
+ cam_out_NETSW   19.729191  38.364034  0.956531
+ cam_out_FLWDS    6.588051   8.529763  0.883237
+ cam_out_PRECSC   3.321645   6.182026      -inf
+ cam_out_PRECC   38.805963  88.711142 -0.926289
+ cam_out_SOLS    10.939536  23.264985  0.929375
+ cam_out_SOLL     13.46133  27.812571  0.915267
+ cam_out_SOLSD    6.158717  12.639114  0.921261
+ cam_out_SOLLD     6.06565  12.624935  0.796327,
+ 'HSR':                       MAE       RMSE         R2
+ variable                                       
+ ptend_t          2.844789   4.824681   0.567524
+ ptend_q0001      4.784426    7.89555       -inf
+ cam_out_NETSW   19.816967  37.773019   0.959136
+ cam_out_FLWDS    6.266635   8.219851   0.904115
+ cam_out_PRECSC   3.511187   6.094945       -inf
+ cam_out_PRECC   42.379195  90.640413 -68.347502
+ cam_out_SOLS    11.314418  23.614026   0.928777
+ cam_out_SOLL    13.602788  27.778323   0.916143
+ cam_out_SOLSD     6.33068   12.39584   0.922621
+ cam_out_SOLLD    6.214983  12.474274   0.797273,
+ 'MLP':                       MAE       RMSE         R2
+ variable                                       
+ ptend_t          2.682765   4.420904   0.588987
+ ptend_q0001      4.494752   7.321821       -inf
+ cam_out_NETSW   13.360914   26.70755   0.982684
+ cam_out_FLWDS    5.224468   6.968642   0.924299
+ cam_out_PRECSC   2.683915    4.73368        NaN
+ cam_out_PRECC   34.333065  72.876357 -38.685671
+ cam_out_SOLS     7.970805  17.403315    0.96138
+ cam_out_SOLL    10.299178  21.954437   0.947612
+ cam_out_SOLSD    4.533141   9.420306   0.955713
+ cam_out_SOLLD    4.806307  10.118702   0.866159,
+ 'RPN':                       MAE       RMSE         R2
+ variable                                       
+ ptend_t          2.685253   4.481825   0.617206
+ ptend_q0001       4.59241   7.518002       -inf
+ cam_out_NETSW   18.878685  33.596484   0.968036
+ cam_out_FLWDS    6.017693   7.913915   0.912173
+ cam_out_PRECSC   3.328444   5.511482       -inf
+ cam_out_PRECC   37.455143  76.578238 -67.943146
+ cam_out_SOLS    10.355971  20.606768   0.942723
+ cam_out_SOLL    12.955313  25.218968   0.928029
+ cam_out_SOLSD    5.845594  10.996472    0.93998
+ cam_out_SOLLD    5.701871  11.247108    0.83663,
+ 'ED':                       MAE       RMSE         R2
+ variable                                       
+ ptend_t          2.864516   4.696491   0.541714
+ ptend_q0001      4.673036   7.642989       -inf
+ cam_out_NETSW   14.967681  28.537017   0.980463
+ cam_out_FLWDS    6.893821   9.070148   0.801872
+ cam_out_PRECSC     3.0459   5.078154       -inf
+ cam_out_PRECC     37.2502  76.681968 -17.908577
+ cam_out_SOLS     8.553734   17.99865   0.959621
+ cam_out_SOLL    10.923703  22.540191   0.945472
+ cam_out_SOLSD    5.074686    9.91694   0.950717
+ cam_out_SOLLD    5.136236  10.416563    0.85718}
+
+
+
+
+
+

Plot 1#

+
+
+
fig, axes = plt.subplots(nrows  = len(data.metrics_names), sharex = True)
+for i in range(len(data.metrics_names)):
+    plot_df_byvar[data.metrics_names[i]].plot.bar(
+        legend = False,
+        ax = axes[i])
+    if data.metrics_names[i] != 'R2':
+        axes[i].set_ylabel('$W/m^2$')
+        axes[i].set_yscale('log')
+    else:
+        axes[i].set_ylim(0,1)
+
+    axes[i].set_title(f'({letters[i]}) {data.metrics_names[i]}')
+axes[i].set_xlabel('Output variable')
+axes[i].set_xticklabels(plot_df_byvar[data.metrics_names[i]].index, \
+    rotation=0, ha='center')
+
+axes[0].legend(columnspacing = .9, 
+               labelspacing = .3,
+               handleheight = .07,
+               handlelength = 1.5,
+               handletextpad = .2,
+               borderpad = .2,
+               ncol = 3,
+               loc = 'upper right')
+fig.set_size_inches(7,8)
+fig.tight_layout()
+
+
+
+
+../_images/54e84fa9f18c91d5f4cb11827da87396e37a4f99c7b678bcc325cb62f7fc5037.png +
+
+
+
+
hmm = data.load_npy_file(load_path = '/ocean/projects/atm200007p/behrens/ED_Behrens_2022/ED_ClimSIM_1_3_pred.npy')
+
+
+
+
+
+
+
hmm.shape
+
+
+
+
+
(1681920, 128)
+
+
+
+
+
+
+
data.target_scoring.shape
+
+
+
+
+
(1681920, 128)
+
+
+
+
+
+
+
import matplotlib
+print(matplotlib.__version__)
+
+
+
+
+
3.5.2
+
+
+
+
+
+
+
fig, ax = plt.subplots()
+ax.plot([1, 2, 3], [4, 5, 6], label='Line 1')
+ax.plot([1, 2, 3], [6, 5, 4], label='Line 2')
+ax.plot([1, 2, 3], [7, 8, 9], label='Line 3')
+ax.legend()
+
+# set the number of columns in the legend to 2
+ax.legend(ncol=2)
+
+# show the plot
+plt.show()
+
+
+
+
+../_images/fbb169e493b54c1fe52192c28f8d9846b852a654b9568f1a9ce1b5dbb0360643.png +
+
+
+
+

Plot 2#

+
+
+

Plot 3#

+
+
+

Plot 4#

+
+
+
+ + + + +
+ + + + + + +
+ + + +
+ +
+
+ Contents +
+
+ +
+ + +
+
+ +
+ +
+ +

+By The ClimSim Maintainers +

+ +
+ +
+ +

+ + © Copyright 2023. +
+ +

+ +
+ +
+ +
+ +
+ +
+ +
+
+ + +
+
+
+ + + + + + + + \ No newline at end of file diff --git a/evaluation/plot_R2_analysis.html b/evaluation/plot_R2_analysis.html new file mode 100644 index 0000000..4ad8289 --- /dev/null +++ b/evaluation/plot_R2_analysis.html @@ -0,0 +1,620 @@ + + + + + + + + + + + + Prediction skill of temperature and specific humidity tendencies — ClimSim + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+
+
+
+ + + Ctrl+K +
+
+ + + +
+
+ +
+ + + +
+ + + + +
+ +
+
+ + +
+
+
+ + +
+
+ + + + +
+ +
+ + + +
+ +
+
+ +
+
+ +
+ +
+ +
+ + +
+ +
+ +
+ + + + + +
+ + +
+ + + + + + + + + + + +
+ +
+ +
+
+ + + +
+

Prediction skill of temperature and specific humidity tendencies

+ +
+ +
+
+ + + + +
+ +
+

Prediction skill of temperature and specific humidity tendencies#

+
+

Importing packages and data_utils.py#

+
+
+
from climsim_utils.data_utils import *
+
+
+
+
+
2023-08-21 16:25:37.644365: I tensorflow/core/platform/cpu_feature_guard.cc:182] This TensorFlow binary is optimized to use available CPU instructions in performance-critical operations.
+To enable the following instructions: AVX2 FMA, in other operations, rebuild TensorFlow with the appropriate compiler flags.
+2023-08-21 16:25:53.040170: W tensorflow/compiler/tf2tensorrt/utils/py_utils.cc:38] TF-TRT Warning: Could not find TensorRT
+
+
+
+
+
+
+

Instantiating class#

+
+
+
grid_path = '/ocean/projects/atm200007p/jlin96/neurips_proj/ClimSim/grid_info/ClimSim_low-res_grid-info.nc'
+norm_path = '/ocean/projects/atm200007p/jlin96/neurips_proj/ClimSim/normalizations/'
+
+grid_info = xr.open_dataset(grid_path)
+input_mean = xr.open_dataset(norm_path + 'inputs/input_mean.nc')
+input_max = xr.open_dataset(norm_path + 'inputs/input_max.nc')
+input_min = xr.open_dataset(norm_path + 'inputs/input_min.nc')
+output_scale = xr.open_dataset(norm_path + 'outputs/output_scale.nc')
+
+scoring_data = data_utils(grid_info = grid_info, 
+                          input_mean = input_mean, 
+                          input_max = input_max, 
+                          input_min = input_min, 
+                          output_scale = output_scale)
+
+
+
+
+
+
+

Loading data and setting pressure grid#

+
+
+
# paths to scoring data
+input_scoring_path = '/ocean/projects/atm200007p/jlin96/neurips_proj/e3sm_train_npy/scoring_input.npy'
+target_scoring_path = '/ocean/projects/atm200007p/jlin96/neurips_proj/e3sm_train_npy/scoring_target.npy'
+
+# paths to model predictions
+cnn_pred_path = '/ocean/projects/atm200007p/jlin96/neurips_proj/figure_ingredients/val_predict_cnn_reshaped_stride6_FINAL.npy'
+ed_pred_path = '/ocean/projects/atm200007p/behrens/ED_Behrens_2022/ED_ClimSIM_1_3_pred.npy'
+hsr_pred_path = '/ocean/projects/atm200007p/shared/neurips_proj/final_metrics/predictions/HSR/hsr_preds_bestcrps.h5'
+mlp_pred_path = '/ocean/projects/atm200007p/jlin96/neurips_proj/figure_ingredients/001_backup_phase-7_retrained_models_step2_lot-147_trial_0027.best.h5.npy'
+rpn_pred_path = '/ocean/projects/atm200007p/jlin96/neurips_proj/figure_ingredients/rpn_pred_v1_stride6.npy'
+cvae_pred_path = '/ocean/projects/atm200007p/shared/neurips_proj/final_metrics/predictions/cVAE/cvae_preds_manual.h5'
+
+# set path to original data
+scoring_data.data_path = '/ocean/projects/atm200007p/jlin96/neurips_proj/e3sm_train/'
+
+# path to target input
+scoring_data.input_scoring = np.load(input_scoring_path)
+
+# path to target output
+scoring_data.target_scoring = np.load(target_scoring_path)
+
+# set pressure weights
+scoring_data.set_pressure_grid(scoring_data.input_scoring)
+
+# set regular expressions for selecting data
+scoring_data.set_regexps(data_split = 'scoring', regexps = ['E3SM-MMF.mli.0008-0[23456789]-*-*.nc', 
+                                                            'E3SM-MMF.mli.0008-1[012]-*-*.nc', 
+                                                            'E3SM-MMF.mli.0009-01-*-*.nc'])
+# set temporal subsampling
+scoring_data.set_stride_sample(data_split = 'scoring', stride_sample = 6)
+
+# create list of files to extract data from
+scoring_data.set_filelist(data_split = 'scoring')
+
+# create pressure grid
+pressure_grid_plotting = scoring_data.get_pressure_grid_plotting(data_split = 'scoring')
+
+# load target output
+scoring_data.target_scoring = np.load(target_scoring_path)
+
+# make a separate copy
+scoring_data_copy = copy.deepcopy(scoring_data)
+
+
+
+
+
100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 4380/4380 [03:04<00:00, 23.72it/s]
+
+
+
+
+
+
+
# load model predictions part 1
+scoring_data.model_names = ['CNN','ED', 'HSR']
+preds = [scoring_data.load_npy_file(load_path = cnn_pred_path),
+         scoring_data.load_npy_file(load_path = ed_pred_path), 
+         scoring_data.load_h5_file(load_path = hsr_pred_path)]
+scoring_data.preds_scoring = dict(zip(scoring_data.model_names, preds))
+
+# load model predictions part 2
+scoring_data_copy.model_names = ['MLP', 'RPN', 'cVAE']
+preds_copy = [scoring_data_copy.load_npy_file(load_path = mlp_pred_path), 
+              scoring_data_copy.load_npy_file(load_path = rpn_pred_path),
+              scoring_data_copy.load_h5_file(load_path = cvae_pred_path)]
+scoring_data_copy.preds_scoring = dict(zip(scoring_data_copy.model_names, preds_copy))
+
+
+
+
+
+
+

Plot pressure latitude R2 analysis#

+
+
+
%config InlineBackend.figure_format = 'retina'
+
+
+
+
+
+
+
scoring_data.plot_r2_analysis(pressure_grid_plotting = pressure_grid_plotting)
+
+
+
+
+
/ocean/projects/atm200007p/jlin96/neurips_proj/ClimSim/climsim_utils/data_utils.py:792: UserWarning: The input coordinates to pcolor are interpreted as cell centers, but are not monotonically increasing or decreasing. This may lead to incorrectly calculated cell edges, in which case, please supply explicit cell edges to pcolor.
+  contour_plot = ax[0,i].pcolor(X, Y, coeff,cmap='Blues', vmin = 0, vmax = 1) # pcolormesh
+/ocean/projects/atm200007p/jlin96/neurips_proj/ClimSim/climsim_utils/data_utils.py:799: RuntimeWarning: divide by zero encountered in divide
+  coeff = 1 - np.sum( (pred_moist_daily_long-test_moist_daily_long)**2, axis=1)/np.sum( (test_moist_daily_long-np.mean(test_moist_daily_long, axis=1)[:,None,:])**2, axis=1)
+/ocean/projects/atm200007p/jlin96/neurips_proj/ClimSim/climsim_utils/data_utils.py:803: UserWarning: The input coordinates to pcolor are interpreted as cell centers, but are not monotonically increasing or decreasing. This may lead to incorrectly calculated cell edges, in which case, please supply explicit cell edges to pcolor.
+  contour_plot = ax[1,i].pcolor(X, Y, coeff,cmap='Blues', vmin = 0, vmax = 1) # pcolormesh
+
+
+../_images/22e0149201f7de3fe4cd14aab3b3af4b18ecd3679df2aac59c8dae9e4301bf5d.png +
<Figure size 640x480 with 0 Axes>
+
+
+
+
+
+
+
scoring_data_copy.plot_r2_analysis(pressure_grid_plotting = pressure_grid_plotting)
+
+
+
+
+
/ocean/projects/atm200007p/jlin96/neurips_proj/ClimSim/climsim_utils/data_utils.py:792: UserWarning: The input coordinates to pcolor are interpreted as cell centers, but are not monotonically increasing or decreasing. This may lead to incorrectly calculated cell edges, in which case, please supply explicit cell edges to pcolor.
+  contour_plot = ax[0,i].pcolor(X, Y, coeff,cmap='Blues', vmin = 0, vmax = 1) # pcolormesh
+/ocean/projects/atm200007p/jlin96/neurips_proj/ClimSim/climsim_utils/data_utils.py:799: RuntimeWarning: divide by zero encountered in divide
+  coeff = 1 - np.sum( (pred_moist_daily_long-test_moist_daily_long)**2, axis=1)/np.sum( (test_moist_daily_long-np.mean(test_moist_daily_long, axis=1)[:,None,:])**2, axis=1)
+/ocean/projects/atm200007p/jlin96/neurips_proj/ClimSim/climsim_utils/data_utils.py:803: UserWarning: The input coordinates to pcolor are interpreted as cell centers, but are not monotonically increasing or decreasing. This may lead to incorrectly calculated cell edges, in which case, please supply explicit cell edges to pcolor.
+  contour_plot = ax[1,i].pcolor(X, Y, coeff,cmap='Blues', vmin = 0, vmax = 1) # pcolormesh
+
+
+../_images/d0a33f5ff37964471db2599b0ed1a252e299637fc3359651f7e46e5e2238a4d8.png +
+
+
+
+ + + + +
+ + + + + + +
+ + + +
+ +
+
+ Contents +
+
+ +
+ + +
+
+ +
+ +
+ +

+By The ClimSim Maintainers +

+ +
+ +
+ +

+ + © Copyright 2023. +
+ +

+ +
+ +
+ +
+ +
+ +
+ +
+
+ + +
+
+
+ + + + + + + + \ No newline at end of file diff --git a/genindex.html b/genindex.html index 5f084f9..45ef9ad 100644 --- a/genindex.html +++ b/genindex.html @@ -97,8 +97,8 @@ class="form-control" name="q" id="search-input" - placeholder="Search..." - aria-label="Search..." + placeholder="Search this book..." + aria-label="Search this book..." autocomplete="off" autocorrect="off" autocapitalize="off" @@ -128,7 +128,7 @@
-