Final update to get_metrics_joint_dom.ipynb.

e4345740 · Michael Langguth · c93a0e2b · e4345740
Commit e4345740 authored 3 years ago by Michael Langguth
--- a/Jupyter_Notebooks/get_metrics_joint_dom.ipynb
+++ b/Jupyter_Notebooks/get_metrics_joint_dom.ipynb
@@ -24,9 +24,12 @@
   "source": [
    "# Evaluation over a smaller (joint) domain\n",
    "\n",
-    "In this Jupyter Notebook, evaluation will be performed on a subregion which is nested into the target region of the evaluated trained video prediction models. The resulting netCDF-files can be used in the meta-postpro\n",
+    "In this Jupyter Notebook, evaluation will be performed on a subregion which is nested into the target region of the evaluated trained video prediction models. The resulting netCDF-files can be used in the meta-postprocessing to perform an evaluation on a joint region such as in Figure 10a of the [manuscript](https://doi.org/10.5194/gmd-2021-430). <br>\n",
-    "The following cells will first merge all forecast files under `indir` into a single netCDF-file.<br>\n",
+    "\n",
-    "Then the data is sliced to the domain defined by `lonlatbox` and all subsequent evaluation is performed on this smaller domain.<br>\n",
+    "### Approach\n",
+    "\n",
+    "In the following cells, all forecast files under `indir` will be merged into a single netCDF-file first.<br>\n",
+    "Then the data gets sliced to the domain defined by `lonlatbox` and all subsequent evaluation is performed on this smaller domain.<br>\n",
    "The evaluation metrics are then saved to a file under `indir` named `evaluation_metrics_<nlon>x<nlat>.nc` where `nlat` and `nlon` denote the number of grid points/pixels in latitude and longitude direction of the smaller domain, respectively. <br>\n",
    "\n",
    "Thus, first let's define the basic parameters:"
@@ -39,9 +42,9 @@
   "metadata": {},
   "outputs": [],
   "source": [
-    "indir = \"/p/project/deepacf/deeprain/video_prediction_shared_folder/results/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/savp/20210901T090059_gong1_savp_cv12/\"\n",
+    "indir = \"<path_to_trained_model/with_larger_target_domain>\"\n",
    "model = \"savp\"\n",
-    "# define domain. [3., 24.3, 40.2, 53.1] corresponds to the smallest domain tested in the GMD paper\n",
+    "# define domain. [3., 24.3, 40.2, 53.1] corresponds to the smallest domain tested in the GMD paper (with 72x44grid points)\n",
    "lonlatbox = [3., 24.3, 40.2, 53.1]"
   ]
  },
@@ -297,9 +300,9 @@
 ],
 "metadata": {
  "kernelspec": {
-   "display_name": "PyDeepLearning-1.0",
+   "display_name": "Python 3",
   "language": "python",
-   "name": "pydeeplearning"
+   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
@@ -311,7 +314,7 @@
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
-   "version": "3.8.5"
+   "version": "3.9.6"
  }
 },
 "nbformat": 4,

 %% Cell type:code id:written-productivity tags:
 ``` python
 import os, sys
 import glob
 sys.path.append("../utils/")
 import xarray as xr
 import numpy as np
 import pandas as pd
 from statistical_evaluation import perform_block_bootstrap_metric, avg_metrics, calculate_cond_quantiles, Scores
 ```
 %% Cell type:markdown id:english-conducting tags:
 # Evaluation over a smaller (joint) domain
-In this Jupyter Notebook, evaluation will be performed on a subregion which is nested into the target region of the evaluated trained video prediction models. The resulting netCDF-files can be used in the meta-postpro
+In this Jupyter Notebook, evaluation will be performed on a subregion which is nested into the target region of the evaluated trained video prediction models. The resulting netCDF-files can be used in the meta-postprocessing to perform an evaluation on a joint region such as in Figure 10a of the [manuscript](https://doi.org/10.5194/gmd-2021-430). <br>
-The following cells will first merge all forecast files under `indir` into a single netCDF-file.<br>
-Then the data is sliced to the domain defined by `lonlatbox` and all subsequent evaluation is performed on this smaller domain.<br>
+### Approach
+In the following cells, all forecast files under `indir` will be merged into a single netCDF-file first.<br>
+Then the data gets sliced to the domain defined by `lonlatbox` and all subsequent evaluation is performed on this smaller domain.<br>
 The evaluation metrics are then saved to a file under `indir` named `evaluation_metrics_<nlon>x<nlat>.nc` where `nlat` and `nlon` denote the number of grid points/pixels in latitude and longitude direction of the smaller domain, respectively. <br>
 Thus, first let's define the basic parameters:
 %% Cell type:code id:metropolitan-mailman tags:
 ``` python
-indir = "/p/project/deepacf/deeprain/video_prediction_shared_folder/results/era5-Y2007-2019M01to12-92x56-3840N0000E-2t_tcc_t_850/savp/20210901T090059_gong1_savp_cv12/"
+indir = "<path_to_trained_model/with_larger_target_domain>"
 model = "savp"
-# define domain. [3., 24.3, 40.2, 53.1] corresponds to the smallest domain tested in the GMD paper
+# define domain. [3., 24.3, 40.2, 53.1] corresponds to the smallest domain tested in the GMD paper (with 72x44grid points)
 lonlatbox = [3., 24.3, 40.2, 53.1]
 ```
 %% Cell type:markdown id:turned-player tags:
 Next, we perform a first merging step. For computational efficiency, we merge max. 1000 files in the first step.<br>
 Since the data is not sorted by the dimension `init_time` when querying along the sample index, we sort it before saving to intermediate files.<br>
 Given that the merging step has already been performed, no further processing is required.<br>
 If this is not the case, we start with the sample indices between 0 and 999:
 %% Cell type:code id:aging-radio tags:
 ``` python
 def get_fname_ndigits(indir, prefix, suffix, n, patt="[0-9]"):
    flist = []
    for i in range(1, n+1):
        fn_search = os.path.join(indir, "{0}{1}{2}".format(prefix, i*patt, suffix))
        flist = flist + glob.glob(fn_search)
    if len(flist) == 0:
        raise FileNotFoundError("Could not find any file under '{0}' with prefix '{1}' and suffix '{2}' containing digits.".format(indir, prefix, suffix))
    return flist
 # get list of files with sample index between 0 and 999.
 vfp_list = get_fname_ndigits(indir, "vfp_date_*sample_ind_", ".nc", 3)
 outfile = os.path.join(indir, "vfp_{0}_forecasts_sample_ind_0_999.nc".format(model))
 if not os.path.isfile(outfile):
    print("File '{0}' does not exist. \n Start reading data with sample index between 0 and 999 from '{1}'...".format(outfile, indir))
    data_all = xr.open_mfdataset(vfp_list, concat_dim="init_time", combine="nested", decode_cf=True).load()
    data_all = data_all.sortby("init_time")
    print("Data loaded successfully. Save merged data to '{0}'.".format(outfile))
    data_all.to_netcdf(outfile, encoding={'init_time':{'units': "seconds since 1900-01-01 00:00:00"}})
 ```
 %% Cell type:markdown id:controversial-picking tags:
 Then, we proceed with the rest.
 %% Cell type:code id:basic-corps tags:
 ``` python
 for i in np.arange(1, 9):
    outfile = os.path.join(indir, "vfp_{0}_forecasts_sample_ind_{1:d}000_{1:d}999.nc".format(model, i))
    if not os.path.isfile(outfile):
        print("File '{0}' does not exist. Start reading data with sample index between {1:d}000 and {1:d}999 from '{2}'...".format(outfile, i, indir))
        data_all = xr.open_mfdataset(os.path.join(indir, "vfp_date_*sample_ind_{0}???.nc".format(i)), concat_dim="init_time", combine="nested", decode_cf=True).load()
        data_all = data_all.sortby("init_time")
        print("Data loaded successfully. Save merged data to '{0}'.".format(outfile))
        data_all.to_netcdf(outfile, encoding={'init_time':{'units': "seconds since 1900-01-01 00:00:00"}})
 ```
 %% Cell type:markdown id:severe-satisfaction tags:
 Still, xarray's `open_mfdataset`-method would not be able to concatenate all data since the `init_time`-dimension is not montonically increasing/decreasing when looping through the files. <br>
 Thus, we have to merge the data manually.
 The merged dataset is then saved to separate datafile for later computation.
 If the data has already been merged, we simply read the data from the corresponding netCDF-file.
 %% Cell type:code id:illegal-headset tags:
 ``` python
 outfile_all = os.path.join(indir, "vfp_{0}_forecasts_all.nc".format(model))
 if not os.path.isfile(outfile_all):
    print("netCDF-file with all forecasts '{0}' does not exist yet. Start merging and sorting all precursor files.".format(outfile))
    all_files = sorted(glob.glob(os.path.join(indir, "vfp_{0}_forecasts_sample_ind_*.nc".format(model))))
    if len(all_files) == 0:
        raise FileNotFoundError("Could not find any precursor files.")
    for i, f in enumerate(all_files):
        print("Processing file '{0}'".format(f))
        tmp = xr.open_dataset(os.path.join(indir, f)).load()
        if i == 0:
            all_fcst = tmp.copy()
        else:
            print("Start merging")
            all_fcst = xr.merge([all_fcst, tmp])
    # sort by init_time-dimension...
    all_fcst = all_fcst.sortby("init_time")
    # ... and save to file
    print("Finally, write all merged and sorted data to '{0}'.".format(outfile_all))
    all_fcst.to_netcdf(outfile_all)
 else:
    all_fcst = xr.open_dataset(outfile_all).load()
 ```
 %% Cell type:markdown id:fresh-favor tags:
 Now, we slice the dataset to the domain of interest (defined by `lonlatbox`).
 %% Cell type:code id:initial-campbell tags:
 ``` python
 all_fcst_sl = all_fcst.sel({"lon": slice(lonlatbox[0], lonlatbox[1]), "lat": slice(lonlatbox[3], lonlatbox[2])})
 print(all_fcst_sl)
 ```
 %% Cell type:markdown id:optimum-drunk tags:
 Next we initialize the function for calculating the MSE and call it to evaluate the ERA5 and persistence forecasts. <br>
 If you require further evaluation metrics, just expand the cell accordingly, e.g. add the following lines <br>
 ```
 ssim_func = Scores("ssim", ["lat", "lon"]).score_func
 ssim_era5_all = ssim_func(data_fcst=era5_fcst[varname_fcst], data_ref=era5_fcst[varname_ref])
 ssim_per_all = (data_fcst=era5_fcst[varname_per], data_ref=era5_fcst[varname_ref])
 ```
 in case you want to evaluate the SSIM as well.
 %% Cell type:code id:sorted-adams tags:
 ``` python
 mse_func = Scores("mse", ["lat", "lon"]).score_func
 varname_ref, varname_fcst, varname_per = "2t_ref", "2t_{0}_fcst".format(model), "2t_persistence_fcst"
 mse_model_all = mse_func(data_fcst=all_fcst_sl[varname_fcst], data_ref=all_fcst_sl[varname_ref])
 mse_per_all = mse_func(data_fcst=all_fcst_sl[varname_per], data_ref=all_fcst_sl[varname_ref])
 ```
 %% Cell type:markdown id:demonstrated-eligibility tags:
 Then, we initialize the data arrays to store the desired evaluation metrics...
 %% Cell type:code id:protective-battle tags:
 ``` python
 fcst_hours = all_fcst_sl["fcst_hour"]
 nhours = len(fcst_hours)
 nboots=1000
 mse_model_fcst = xr.DataArray(np.empty(nhours, dtype=object), coords={"fcst_hour": fcst_hours}, dims=["fcst_hour"])
 mse_model_fcst_boot = xr.DataArray(np.empty((nhours, nboots), dtype=object),
                                  coords={"fcst_hour": fcst_hours, "iboot": np.arange(nboots)},
                                  dims=["fcst_hour", "iboot"])
 mse_per_fcst = xr.DataArray(np.empty(nhours, dtype=object), coords={"fcst_hour": fcst_hours}, dims=["fcst_hour"])
 mse_per_fcst_boot = xr.DataArray(np.empty((nhours, nboots), dtype=object),
                                 coords={"fcst_hour": fcst_hours, "iboot": np.arange(nboots)},
                                 dims=["fcst_hour", "iboot"])
 ```
 %% Cell type:markdown id:norman-swedish tags:
 ... and populate them by looping over all forecast hours.
 %% Cell type:code id:extra-bible tags:
 ``` python
 for i, fh in enumerate(fcst_hours):
    mse_model_curr = mse_model_all.sel(fcst_hour=fh)
    mse_per_curr = mse_per_all.sel(fcst_hour=fh)
    mse_model_fcst[fh-1], mse_per_fcst[fh-1] = mse_model_curr.mean(), mse_per_curr.mean()
    mse_model_fcst_boot[i, :] = perform_block_bootstrap_metric(mse_model_curr, "init_time", 24*7)
    mse_per_fcst_boot[i, :] = perform_block_bootstrap_metric(mse_per_curr, "init_time", 24*7)
 ```
 %% Cell type:markdown id:limiting-corpus tags:
 Finally, we put the data arrays into a joint dataset and save the results into the netCDF-file.
 %% Cell type:code id:reported-traffic tags:
 ``` python
 # create Dataset and save to netCDF-file
 ds_mse = xr.Dataset({"2t_{0}_mse_avg".format(model): mse_model_fcst, "2t_{0}_mse_bootstrapped".format(model): mse_model_fcst_boot,
                     "2t_persistence_mse_avg": mse_per_fcst, "2t_persistence_mse_bootstrapped": mse_per_fcst_boot})
 outfile = os.path.join(indir, "evaluation_metrics_{0:d}x{1:d}.nc".format(len(all_fcst_sl["lon"]), len(all_fcst_sl["lat"])))
 print("Save evaluation metrics to '{0}'".format(outfile))
 print(ds_mse)
 ds_mse.to_netcdf(outfile)
 ```
 %% Cell type:markdown id:correct-bulgarian tags:
 ## Done!