Merge branch 'master' into housekeeping

AusClimateService · Dec 6, 2021 · 1662221 · 1662221
2 parents d4d4672 + cc1ac73
commit 1662221
Show file tree

Hide file tree

Showing 13 changed files with 2,917 additions and 645 deletions.
diff --git a/cmdline_scripts/bias_correct.py b/cmdline_scripts/bias_correct.py
@@ -22,21 +22,19 @@
 def _main(args):
     """Run the command line program."""
 
-    ds_fcst = fileio.open_file(args.fcst_file)
-    da_fcst = ds_fcst[args.var]
-    # init_dates = time_utils.cftime_to_str(da_fcst["init_date"])
-    # n_lead_steps = int(da_fcst["lead_time"].values.max()) + 1
-
     ds_obs = fileio.open_file(args.obs_file, variables=[args.var])
     da_obs = ds_obs[args.var]
 
+    ds_fcst = fileio.open_file(args.fcst_file, variables=[args.var])
+    da_fcst = ds_fcst[args.var]
+    init_dates = time_utils.cftime_to_str(da_fcst["init_date"])
+    n_lead_steps = int(da_fcst["lead_time"].values.max()) + 1
+
     bias = bias_correction.get_bias(
-        da_fcst, da_obs, args.method, time_period=args.base_period, monthly=args.monthly
-    )
-    da_fcst_bc = bias_correction.remove_bias(
-        da_fcst, bias, args.method, monthly=args.monthly
+        da_fcst, da_obs, args.method, time_period=args.base_period
     )
-
+    da_fcst_bc = bias_correction.remove_bias(da_fcst, bias, args.method)
+
     ds_fcst_bc = da_fcst_bc.to_dataset()
     infile_logs = {
         args.fcst_file: ds_fcst.attrs["history"],
@@ -73,12 +71,6 @@ def _main(args):
         nargs=2,
         help="Start and end date for baseline (YYYY-MM-DD format)",
     )
-    parser.add_argument(
-        "--monthly",
-        action="store_true",
-        default=False,
-        help="Monthly climatology / bias removal [default=False]",
-    )
     parser.add_argument(
         "--output_chunks",
         type=str,

diff --git a/projects/seatac/Makefile b/projects/seatac/Makefile
@@ -1,13 +1,54 @@
-.PHONY: all plot-histogram plot-reanalysis-hot-day plot-model-hot-day plot-sample-size-dist plot-likelihoods plot-return-periods plot-annual-max plot-distribution clean help
+.PHONY: all tasmax-obs tasmax-forecast tasmax-bias-correction calc-txx-obs calc-txx-forecast calc-txx-forecast-bias-corrected plot-histogram plot-reanalysis-hot-day plot-model-hot-day plot-sample-size-dist plot-likelihoods plot-return-periods plot-annual-max plot-distribution clean help
 
 include ${CONFIG}
 
 PYTHON=/g/data/xv83/dbi599/miniconda3/envs/unseen/bin/python
+SCRIPT_DIR=/home/599/dbi599/unseen/cmdline_scripts
 PLOT_PARAMS=${CONFIG_DIR}/plotparams_publication.yml
 
+## tasmax-obs : preparation of observed tasmax data
+tasmax-obs : ${OBS_TASMAX_FILE}
+${OBS_TASMAX_FILE} : ${OBS_DATA} ${OBS_METADATA}
+	${PYTHON} ${SCRIPT_DIR}/preprocess.py $< obs $@ --metadata_file $(word 2,$^) --variables tasmax --units ${UNITS} --no_leap_day --input_freq D 
+
+## tasmax-forecast : preparation of CAFE tasmax data
+tasmax-forecast : ${FCST_TASMAX_FILE}
+${FCST_TASMAX_FILE} : ${FCST_METADATA}
+	${PYTHON} ${SCRIPT_DIR}/preprocess.py ${FCST_DATA} forecast $@ --metadata_file $< --variables tasmax --units ${UNITS} --no_leap_day --input_freq D --spatial_coords ${LAT} ${LON} --output_chunks lead_time=50 --dask_config ${DASK_CONFIG}
+
+## tasmax-bias-correction : bias correct tasmax data using observations
+tasmax-bias-correction : ${FCST_TASMAX_BIAS_CORRECTED_FILE}
+${FCST_TASMAX_BIAS_CORRECTED_FILE} : ${FCST_TASMAX_FILE} ${OBS_TASMAX_FILE}
+	${PYTHON} ${SCRIPT_DIR}/bias_correct.py $< $(word 2,$^) tasmax ${BIAS_METHOD} $@ --base_period ${BASE_PERIOD}
+
+## calc-txx-obs : calculate annual daily maximum temperature from observational data
+calc-txx-obs : ${OBS_TXX_FILE}
+${OBS_TXX_FILE} : ${OBS_TASMAX_FILE}
+	${PYTHON} ${SCRIPT_DIR}/preprocess.py $< obs $@ --variables tasmax --time_freq ${TIME_FREQ} --time_agg ${TIME_AGG} --input_freq D 
+
+## calc-txx-forecast : calculate annual daily maximum temperature from forecast data
+calc-txx-forecast : ${FCST_TXX_FILE}
+${FCST_TXX_FILE} : ${FCST_TASMAX_FILE}
+	${PYTHON} ${SCRIPT_DIR}/preprocess.py $< obs $@ --variables tasmax --time_freq ${TIME_FREQ} --time_agg ${TIME_AGG} --complete_time_agg_periods --input_freq D 
+
+## calc-txx-forecast-bias-corrected : calculate annual daily maximum temperature from bias corrected forecast data
+calc-txx-forecast-bias-corrected : ${FCST_TXX_BIAS_CORRECTED_FILE}
+${FCST_TXX_BIAS_CORRECTED_FILE} : ${FCST_TASMAX_BIAS_CORRECTED_FILE}
+	${PYTHON} ${SCRIPT_DIR}/preprocess.py $< obs $@ --variables tasmax --time_freq ${TIME_FREQ} --time_agg ${TIME_AGG} --complete_time_agg_periods --input_freq D 
+
+## similarity-test : similarity test between observations and bias corrected forecast
+#similarity-test : ${SIMILARITY_FILE}
+#${SIMILARITY_FILE} : ${FCST_BIAS_FILE} ${OBS_PROCESSED_FILE}
+#	${PYTHON} ${SCRIPT_DIR}/similarity_test.py $< $(word 2,$^) ${VAR} $@ --reference_time_period ${BASE_PERIOD}
+
+## independence-test : independence test for different lead times
+#independence-test : ${INDEPENDENCE_PLOT}
+#${INDEPENDENCE_PLOT} : ${FCST_BIAS_FILE}
+#	${PYTHON} ${SCRIPT_DIR}/independence_test.py $< ${VAR} $@ ${INDEPENDENCE_OPTIONS}
+
 ## plot-histogram : plot TXx histogram
 plot-histogram : ${TXX_HISTOGRAM_PLOT}
-${TXX_HISTOGRAM_PLOT} : ${OBS_PROCESSED_FILE} ${FCST_ENSEMBLE_FILE} ${FCST_BIAS_FILE}
+${TXX_HISTOGRAM_PLOT} : ${OBS_TXX_FILE} ${FCST_TXX_FILE} ${FCST_TXX_BIAS_CORRECTED_FILE}
 	${PYTHON} plot_seatac_TXx_histogram.py $< $(word 2,$^) $(word 3,$^) $@  --plotparams ${PLOT_PARAMS}
 
 ## plot-reanalysis-hot-day : plot reanalysis hottest day
@@ -18,31 +59,31 @@ ${REANALYSIS_HOT_DAY_PLOT} : ${REANALYSIS_HGT_FILE} ${REANALYSIS_TAS_FILE}
 ## plot-model-hot-day : plot model hottest day
 plot-model-hot-day : ${MODEL_HOT_DAY_PLOT}
 ${MODEL_HOT_DAY_PLOT} : ${FCST_HOT_DAY_DATA} ${FCST_METADATA}
-	${PYTHON} plot_model_hottest_day.py $< $(word 2,$^) ${LAT} ${LON} $@ --plotparams ${PLOT_PARAMS}
+	${PYTHON} plot_model_hottest_day.py $< $(word 2,$^) ${LAT} ${LON} ${FCST_HOT_DAY_YEAR} $@ --plotparams ${PLOT_PARAMS}
 
 ## plot-sample-size-dist : plot TXx sample size distribution
 plot-sample-size-dist : ${TXX_SAMPLE_PLOT}
-${TXX_SAMPLE_PLOT} : ${FCST_BIAS_FILE}
+${TXX_SAMPLE_PLOT} : ${FCST_TXX_BIAS_CORRECTED_FILE}
 	${PYTHON} plot_TXx_sample_size_distribution.py $< $@ --plotparams ${PLOT_PARAMS}
 
 ## plot-likelihoods : plot TXx likelihoods
 plot-likelihoods : ${TXX_LIKELIHOOD_PLOT}
-${TXX_LIKELIHOOD_PLOT} : ${FCST_BIAS_FILE}
+${TXX_LIKELIHOOD_PLOT} : ${FCST_TXX_BIAS_CORRECTED_FILE}
 	${PYTHON} plot_TXx_likelihoods.py $< $@ --plotparams ${PLOT_PARAMS}
 
 ## plot-return-periods : plot TXx return periods
 plot-return-periods : ${TXX_RETURN_PERIODS_PLOT}
-${TXX_RETURN_PERIODS_PLOT} : ${FCST_BIAS_FILE}
+${TXX_RETURN_PERIODS_PLOT} : ${FCST_TXX_BIAS_CORRECTED_FILE}
 	${PYTHON} plot_TXx_return_periods.py $< $@ --plotparams ${PLOT_PARAMS}
 
 ## plot-annual-max : plot maximum TXx by year
 plot-annual-max : ${TXX_ANNUAL_MAX_PLOT}
-${TXX_ANNUAL_MAX_PLOT} : ${FCST_BIAS_FILE}
+${TXX_ANNUAL_MAX_PLOT} : ${FCST_TXX_BIAS_CORRECTED_FILE}
 	${PYTHON} plot_maxTXx_by_year.py $< $@ --plotparams ${PLOT_PARAMS}
 
 ## plot-distribution : plot TXx distribution by year
 plot-distribution : ${TXX_ANNUAL_DIST_PLOT}
-${TXX_ANNUAL_DIST_PLOT} : ${FCST_BIAS_FILE}
+${TXX_ANNUAL_DIST_PLOT} : ${FCST_TXX_BIAS_CORRECTED_FILE}
 	${PYTHON} plot_TXx_distribution.py $< $@ --plotparams ${PLOT_PARAMS}
 
 ## clean : remove all generated files

diff --git a/projects/seatac/README.md b/projects/seatac/README.md
@@ -4,13 +4,7 @@ Risbey JS, Irving DB, Squire DT, Matear RJ, Monselesan DP, Richardson D & Tozer
 On the role of weather and sampling in assessing a record-breaking heat extreme.
 *Environmental Research Letters*.
 
-The biased corrected forecast ensemble file was created by running the `Makefile` in the root directory of this repository:
-```bash
-$ make bias-correction CONFIG=projects/seatac/config_seatac.mk
-```
-
-The `projects/seatac/Makefile` then defines the rules to use that ensemble file to
-make most of the figures in the paper.
+The `Makefile` defines the rules to make the figures in the paper.
 - Figure 2: `make plot-historgram CONFIG=config_seatac.mk`
 - Figure 4: `make plot-reanalysis-hot-day CONFIG=config_seatac.mk`
 - Figure 5: `make plot-model-hot-day CONFIG=config_seatac.mk`

diff --git a/projects/seatac/config_seatac.mk b/projects/seatac/config_seatac.mk
@@ -1,37 +1,36 @@
 # Configuration for SeaTac airport annual maximum temperture analysis
 
-VAR=tasmax
-UNITS=${VAR}=C
+UNITS=tasmax=C
 LAT=47.45
 LON=237.69
 REGION_NAME=seatac
 TIME_FREQ=A-DEC
 TIME_AGG=max
 BIAS_METHOD=additive
-BASE_PERIOD=1995-01-01 2021-12-31
+BASE_PERIOD=2004-01-01 2020-12-31
 
 CONFIG_DIR=/home/599/dbi599/unseen/config
 SEATAC_DIR=/g/data/xv83/dbi599/seatac
 
 DASK_CONFIG=${CONFIG_DIR}/dask_local.yml
 
-GENERAL_IO_OPTIONS=--variables ${VAR} --time_freq ${TIME_FREQ} --time_agg ${TIME_AGG} --units ${UNITS}
-FCST_IO_OPTIONS=${GENERAL_IO_OPTIONS} --spatial_coords ${LAT} ${LON} --complete_time_agg_periods
-OBS_IO_OPTIONS=${GENERAL_IO_OPTIONS} --input_freq D 
-
 FCST_DATA_1990S := $(sort $(wildcard /g/data/xv83/dcfp/CAFE-f6/c5-d60-pX-f6-199[5,6,7,8,9]*/atmos_isobaric_daily.zarr.zip))
 FCST_DATA_2000S := $(sort $(wildcard /g/data/xv83/dcfp/CAFE-f6/c5-d60-pX-f6-2*/atmos_isobaric_daily.zarr.zip))
 FCST_DATA := ${FCST_DATA_1990S} ${FCST_DATA_2000S}
-FCST_HOT_DAY_DATA = /g/data/xv83/dcfp/CAFE-f6/c5-d60-pX-f6-20110501/atmos_isobaric_daily.zarr.zip
+FCST_HOT_DAY_DATA = /g/data/xv83/dcfp/CAFE-f6/c5-d60-pX-f6-20021101/atmos_isobaric_daily.zarr.zip
+FCST_HOT_DAY_YEAR = 2003
 FCST_METADATA=${CONFIG_DIR}/dataset_cafe_daily.yml
-FCST_ENSEMBLE_FILE=${SEATAC_DIR}/data/${VAR}_cafe-c5-d60-pX-f6_19950501-20201101_${TIME_FREQ}-${TIME_AGG}_${REGION_NAME}.zarr.zip
 FCST_CONFIG=${CONFIG_DIR}/dataset_cafe_daily.yml
+FCST_TASMAX_FILE=${SEATAC_DIR}/data/tasmax_cafe-c5-d60-pX-f6_19950501-20201101_${REGION_NAME}.zarr.zip
+FCST_TXX_FILE=${SEATAC_DIR}/data/tasmax_cafe-c5-d60-pX-f6_19950501-20201101_${TIME_FREQ}-${TIME_AGG}_${REGION_NAME}.zarr.zip
 
 OBS_DATA=${SEATAC_DIR}/data/tasmax_GHCNv2_SeaTac_daily_1948-2021.nc
 OBS_METADATA=${CONFIG_DIR}/dataset_knmi_daily.yml
-OBS_PROCESSED_FILE=${SEATAC_DIR}/data/${VAR}_GHCNv2_1948-2021_${TIME_FREQ}-${TIME_AGG}_${REGION_NAME}.zarr.zip
+OBS_TASMAX_FILE=${SEATAC_DIR}/data/tasmax_GHCNv2_1948-2021_${REGION_NAME}.zarr.zip
+OBS_TXX_FILE=${SEATAC_DIR}/data/tasmax_GHCNv2_1948-2021_${TIME_FREQ}-${TIME_AGG}_${REGION_NAME}.zarr.zip
 
-FCST_BIAS_FILE=${SEATAC_DIR}/data/${VAR}_cafe-c5-d60-pX-f6_19950501-20201101_${TIME_FREQ}-${TIME_AGG}_${REGION_NAME}_bias-corrected-GHCNv2-${BIAS_METHOD}.zarr.zip
+FCST_TASMAX_BIAS_CORRECTED_FILE=${SEATAC_DIR}/data/tasmax_cafe-c5-d60-pX-f6_19950501-20201101_${REGION_NAME}_bias-corrected-GHCNv2-${BIAS_METHOD}.zarr.zip
+FCST_TXX_BIAS_CORRECTED_FILE=${SEATAC_DIR}/data/tasmax_cafe-c5-d60-pX-f6_19950501-20201101_${TIME_FREQ}-${TIME_AGG}_${REGION_NAME}_bias-corrected-GHCNv2-${BIAS_METHOD}.zarr.zip
 
 REANALYSIS_HGT_FILE=${SEATAC_DIR}/data/h500_ERA5_hourly_2021-06-28.grib
 REANALYSIS_TAS_FILE=${SEATAC_DIR}/data/tas_ERA5_hourly_2021-06-28.grib