Produce NetCDF file with 3-day wet hour sums

SEPTREFHUM.py

@@ -6,6 +6,10 @@ from datetime import datetime, timedelta
 # TODO: Figure out where files will be collected from and stored to.
 infile_path  = "in/"
 outfile_path = "out/"
+
+# Remove all out  files
+subprocess.run("rm %s*.nc" % outfile_path, shell=True)
+
 weatherdata_files = glob.glob("%sanalysis-*.nc" % infile_path)
 
 # Iterate the set of hourly weather data files
@@ -17,15 +21,58 @@ for file in weatherdata_files:
     # Produce daily files with WH_SUM, which is the number of "Wet hours" (WH) for a given day
     # WH is defined as RR > 0.2 || UM > 88.0
     subprocess.run(
-        'cdo -setdate,%s -settime,22:00:00 -chname,WH,WH_SUM -timsum -selname,WH -aexpr,"WH = RR > 0.2 || UM > 88.0 ? 1 : 0;" %s %swh_%s.nc' % (wh_sum_date.strftime("%Y-%m-%d"), file, outfile_path, file_date),
+        'cdo -setdate,%s -settime,22:00:00 -chname,WH,WH_DAYSUM -timsum -selname,WH -aexpr,"WH = RR > 0.2 || UM > 88.0 ? 1 : 0;" %s %swh_%s.nc' % (wh_sum_date.strftime("%Y-%m-%d"), file, outfile_path, file_date),
         shell=True
         )
 
+
 # Concatenate daily files > one file with daily values
 subprocess.run(
-        'cdo -O mergetime %swh_*.nc %swh.nc' % (outfile_path, outfile_path),
+        'cdo -O mergetime %swh_*.nc %swh_daysum.nc' % (outfile_path, outfile_path),
         shell=True
         )
 
-# TODO: Add sum of WH_SUM[yesterday] + WH_SUM[today] + WH_SUM[tomorrow]
+# Add sum of WH_SUM[yesterday] + WH_SUM[today] + WH_SUM[tomorrow]
+# timselsum skips every 3 steps when summing 3 timestemps, so we must
+# create three different files and then merge them
+subprocess.run('cdo timselsum,3,0 out/wh_daysum.nc out/wh_3daysum_tmp_0.nc', shell=True)
+subprocess.run('cdo timselsum,3,1 out/wh_daysum.nc out/wh_3daysum_tmp_1.nc', shell=True)
+subprocess.run('cdo timselsum,3,2 out/wh_daysum.nc out/wh_3daysum_tmp_2.nc', shell=True)
+
+subprocess.run('cdo -mergetime out/wh_3daysum_tmp_*.nc out/wh_3daysum_tmp_merged.nc', shell=True)
+
+# the last timesteps are most likely wrong, due to lack of "tomorrows" when performing timselsum
+# To remove the last ones:
+#
+# The variable time_bnds (which has two dimensions: time and bnds. bnds is size 2) contains the 
+# time (in seconds) between the first and last timesteps used for summing up. e.g like this:
+# $ ncdump -v time_bnds out/wh_3daysum.nc
+# time_bnds =
+# 1696111200, 1696284000,
+# 1696197600, 1696370400,
+# 1696284000, 1696456800,
+# 1696370400, 1696543200,
+# 1696456800, 1696629600,
+# 1696543200, 1696716000,
+# 1696629600, 1696716000,
+# 1696716000, 1696716000;
+#}
+# The difference [1] - [0] should be 172800 seconds = 47 hours
+# Timesteps with [1] - [0] != 172800 should be discarded
+# Using netCDF4 to accomplish this
+
+import netCDF4 as nc
+import numpy as np
+
+
+wh_3daysum = nc.Dataset('out/wh_3daysum_tmp_merged.nc', 'r')
+time_bnds = wh_3daysum.variables["time_bnds"][:] 
+# Assuming that wrong time bounds only exist at the end of the time series, this works
+number_of_timesteps_to_remove = 0
+for time_bnd in time_bnds:
+    if time_bnd[1]-time_bnd[0] != 172800.0:
+        number_of_timesteps_to_remove = number_of_timesteps_to_remove + 1
+wh_3daysum.close()
+number_of_timesteps_to_keep = len(time_bnds) - number_of_timesteps_to_remove
+subprocess.run('cdo -seltimestep,1/%s out/wh_3daysum_tmp_merged.nc out/wh_3daysum.nc' % number_of_timesteps_to_keep, shell=True)