Merge pull request #75 from ISI-MIP/develop

Develop

Author: SimonTreu

README.md

@@ -83,67 +83,50 @@ You may optionally
 
 ## Comments for each variable
 
-#### tas
-data checked
-Works using Normal distribution
+#### daily mean temperature (tas)
+Two cells fail in complete dataset.
+67418 of 6420 cells.
 
-#### rlds
-data checked
-Works using Normal distribution
-Needs a restart to finish some hanging runs
+#### tasskew
+Work in progress. See #59
 
-#### psl / ps
-data checked
-Works using Normal distribution
+#### tasrange
+Calculation alsmost complete on full dataset.
+Some cells do not detrend as expected. Need assessment.
+See #60.
 
-#### rsds
-Deviationg approach from Lange et al. 2019, using Normal distribution
-This is because the yearly cycle is handled inherently here, so no need for specific treatment.
-FIXME: produces unrealistic incoming radiation below zero. Needs a different approach
+#### precipiation (pr)
+Calculatio complete.
+67339 of 67420 work.
 
-#### hurs (relative humidity)
-data checked
-With Beta distribution, working
-Needs to be rerun so some holes are filled.
+#### sea level pressure (ps)
+Calculation complete on full dataset.
+See #65
 
-GSWP: needs preprocessing to rename from rhs to hurs, and mask invalid values below zero:
+#### wind
+Calculation complete on full dataset.
+Minor issues on the coast of the Arabic Peninsula.
+See #66
 
-```
-ncrename -O -v rhs,hurs fname1.nc fname2.nc
+#### longwave radiation (rlds)
+Calculation complete on full dataset.
+Need check of trend removal.
 
-cdo setrtomiss,-1e20,0 fname2.nc fname3.nc
-```
+#### shortwave radiation (rsds)
 
-#### tasskew
-data checked
-Works using Beta distribution
+#### relative humidity (hurs)
 
-#### prsnratio
-Beta distribution
-Snow included in GSWP3
 
-#### tasrange
-With Rice distribution
-ADVI introduces strong positive trend.
-Possible issue: use real mu, not nu for quantile mapping.
 
-#### tasmin
-Constructed from tas, tasskew and tasrange
-To do in postprocessing
+## Comments for datasets
 
-#### tasmax
-Constructed from tas, tasskew and tasrange
-To do in postprocessing
+GSWP: needs preprocessing to rename from rhs to hurs, and mask invalid values below zero:
 
-#### pr
-Gamma distribution
-Does not remove all regional trends with NUTS.
-Fails with ADVI.
-Low latitudes are particularly most difficult.
+```
+ncrename -O -v rhs,hurs fname1.nc fname2.nc
 
-#### wind
-Works using Weibull distribution
-FIXME: does not seem to detrend. Seems we rather chose the parameter that adjusted the variability range
+cdo setrtomiss,-1e20,0 fname2.nc fname3.nc
+```
 
 
 ## Credits

icounter/const.py

@@ -23,11 +23,13 @@
     "prsnratio": (0.0, 1.0),
     "hurs": (0, 100.0),
     "ps": (0, None),
-    "rsds": (0, None),
+    "rsds": (0, 501),
     "rlds": (0, None),
     "wind": (0.0, None),
 }
 
+# for rsds, the global maximum is 500
+
 
 def check_bounds(data, variable):
 
@@ -128,7 +130,7 @@ def rescale_and_offset_precip(scaled_data, datamin, scale):
     "tas": [scale_to_unity, rescale_to_original],
     "ps": [scale_to_unity, rescale_to_original],
     "rlds": [scale_to_unity, rescale_to_original],
-    "rsds": [scale_to_unity, rescale_to_original],
+    "rsds": [mask_and_scale_by_bounds, refill_and_rescale],
     "wind": [scale_and_mask, refill_and_rescale],
     "hurs": [mask_and_scale_by_bounds, refill_and_rescale],
     "prsnratio": [mask_and_scale_by_bounds, refill_and_rescale],

icounter/datahandler.py

@@ -15,7 +15,7 @@ def create_output_dirs(output_dir):
         (output_dir / d).mkdir(parents=True, exist_ok=True)
 
 
-def make_cell_output_dir(output_dir, sub_dir, lat, lon, variable=None):
+def make_cell_output_dir(output_dir, sub_dir, lat, lon, variable):
 
     """ params: output_dir: a pathlib object """
 
@@ -29,7 +29,7 @@ def make_cell_output_dir(output_dir, sub_dir, lat, lon, variable=None):
         return lat_sub_dir
 
 
-def get_valid_subset(df, subset, seed):
+def get_subset(df, subset, seed):
 
     orig_len = len(df)
     if subset > 1:
@@ -38,23 +38,10 @@ def get_valid_subset(df, subset, seed):
         df = df.loc[np.sort(subselect), :].copy()
 
     df.replace([np.inf, -np.inf], np.nan, inplace=True)
-    df_valid = df.dropna(axis=0, how="any")
 
-    print(len(df_valid), "data points used from originally", orig_len, "datapoints.")
+    print(len(df), "data points used from originally", orig_len, "datapoints.")
 
-    return df_valid
-
-
-# def get_valid_index(df, subset, seed):
-
-#     orig_len = len(df)
-#     if subset > 1:
-#         np.random.seed(seed)
-#         subselect = np.random.choice(orig_len, np.int(orig_len/subset), replace=False)
-#         df = df.loc[np.sort(subselect), :].copy()
-
-#     df.replace([np.inf, -np.inf], np.nan, inplace=True)
-#     return df.dropna(axis=0, how="any").index
+    return df
 
 
 def create_dataframe(nct_array, units, data_to_detrend, gmt, variable):
@@ -101,40 +88,26 @@ def create_dataframe(nct_array, units, data_to_detrend, gmt, variable):
     return tdf, datamin, scale
 
 
-def create_ref_df(df, trace_for_qm, ref_period, scale_variability, is_precip=False):
+def create_ref_df(df, trace_for_qm, ref_period, params):
 
     df_params = pd.DataFrame(index=df.index)
 
-    # print(trace_for_qm["mu"])
-
-    df_params.loc[:, "mu"] = trace_for_qm["mu"].mean(axis=0)
-    df_params.loc[:, "sigma"] = trace_for_qm["sigma"].mean(axis=0)
-    if is_precip:
-        df_params.loc[:, "pbern"] = trace_for_qm["pbern"].mean(axis=0)
+    for p in params:
+        df_params.loc[:, p] = trace_for_qm[p].mean(axis=0)
 
     df_params.index = df["ds"]
 
     df_params_ref = df_params.loc[ref_period[0] : ref_period[1]]
     # mean over all years for each day
     df_params_ref = df_params_ref.groupby(df_params_ref.index.dayofyear).mean()
 
-    # case of not scaling variability
-    df_params.loc[:, "sigma_ref"] = df_params["sigma"]
     # write the average values for the reference period to each day of the
     # whole timeseries
     for day in df_params_ref.index:
-        df_params.loc[df_params.index.dayofyear == day, "mu_ref"] = df_params_ref.loc[
-            day, "mu"
-        ]
-        if is_precip:
+        for p in params:
             df_params.loc[
-                df_params.index.dayofyear == day, "pbern_ref"
-            ] = df_params_ref.loc[day, "pbern"]
-        # case of scaling sigma
-        if scale_variability:
-            df_params.loc[
-                df_params.index.dayofyear == day, "sigma_ref"
-            ] = df_params_ref.loc[day, "sigma"]
+                df_params.index.dayofyear == day, p + "_ref"
+            ] = df_params_ref.loc[day, p]
 
     return df_params
 
@@ -161,64 +134,36 @@ def get_source_timeseries(data_dir, dataset, qualifier, variable, lat, lon):
     obs_data.close()
     return df
 
+def get_cell_filename(outdir_for_cell, lat, lon, settings):
 
-def save_to_disk(df_with_cfact, settings, lat, lon, dformat=".h5"):
-
-    outdir_for_cell = make_cell_output_dir(
-        settings.output_dir, "timeseries", lat, lon, settings.variable
+    return outdir_for_cell / (
+        "ts_" + settings.dataset + "_lat" + str(lat) + "_lon" + str(lon) + settings.storage_format
     )
 
-    fname = outdir_for_cell / (
-        "ts_" + settings.dataset + "_lat" + str(lat) + "_lon" + str(lon) + dformat
-    )
+def test_if_data_valid_exists(fname):
 
-    if dformat == ".csv":
-        df_with_cfact.to_csv(fname)
-    elif dformat == ".h5":
-        df_with_cfact.to_hdf(fname, "lat_" + str(lat) + "_lon_" + str(lon), mode="w")
+    if ".h5" in str(fname):
+        pd.read_hdf(fname)
+    elif ".csv" in str(fname):
+        pd.read_csv(fname)
     else:
-        raise NotImplementedError("choose storage format .h5 or csv.")
+        raise ValueError
 
-    print("Saved timeseries to ", fname)
+def save_to_disk(df_with_cfact, fname, lat, lon, storage_format):
 
+    # outdir_for_cell = make_cell_output_dir(
+    #     settings.output_dir, "timeseries", lat, lon, settings.variable
+    # )
 
-def read_from_disk(data_path):
+    # fname = outdir_for_cell / (
+    #     "ts_" + settings.dataset + "_lat" + str(lat) + "_lon" + str(lon) + dformat
+    # )
 
-    if data_path.split(".")[-1] == "h5":
-        df = pd.read_hdf(data_path)
-    elif data_path.split(".")[-1] == "csv":
-        df = pd.read_csv(data_path, index_col=0)
+    if storage_format == ".csv":
+        df_with_cfact.to_csv(fname)
+    elif storage_format == ".h5":
+        df_with_cfact.to_hdf(fname, "lat_" + str(lat) + "_lon_" + str(lon), mode="w")
     else:
         raise NotImplementedError("choose storage format .h5 or csv.")
 
-    return df
-
-
-def form_global_nc(ds, time, lat, lon, vnames, torigin):
-
-    ds.createDimension("time", None)
-    ds.createDimension("lat", lat.shape[0])
-    ds.createDimension("lon", lon.shape[0])
-
-    times = ds.createVariable("time", "f8", ("time",))
-    longitudes = ds.createVariable("lon", "f8", ("lon",))
-    latitudes = ds.createVariable("lat", "f8", ("lat",))
-    for var in vnames:
-        data = ds.createVariable(
-            var,
-            "f4",
-            ("time", "lat", "lon"),
-            chunksizes=(time.shape[0], 1, 1),
-            fill_value=1e20,
-        )
-    times.units = torigin
-    latitudes.units = "degree_north"
-    latitudes.long_name = "latitude"
-    latitudes.standard_name = "latitude"
-    longitudes.units = "degree_east"
-    longitudes.long_name = "longitude"
-    longitudes.standard_name = "longitude"
-    # FIXME: make flexible or implement loading from source data
-    latitudes[:] = lat
-    longitudes[:] = lon
-    times[:] = time
+    print("Saved timeseries to ", fname)