Issue #64: run with subdaily grids

.gitignore

@@ -18,7 +18,10 @@ lisflood_lisvap.egg-info/
 !tests/data/tests_cordex.xml
 !/tests/data/tests_glofas.xml
 !tests/data/tests_efas_1arcmin.xml
+!tests/data/tests_efas_1arcmin_6hourly.xml
+!tests/data/tests_efas_1arcmin_hourly.xml
 !tests/data/tests_efas_1arcmin_yearly.xml
+!tests/data/tests_efas_1arcmin_yearly_output.xml
 !tests/data/tests_efas_1arcmin_360days_calendar.xml
 
 # unignoring reference and input data for tests

settings_tpl.xml

@@ -35,7 +35,7 @@ You can use $(SettingsDir) or $(SettingsPath) to refer directory containing the
 
             <textvar name="DtSec" value="86400">
                 <comment>
-                    time step [seconds] ALWAYS USE 86400!!
+                    time step [seconds] USE 86400 for daily, 43200 for 12hourly, 21600 for 6hourly and 3600 for hourly!!
                 </comment>
             </textvar>
 

src/lisvap/__init__.py

@@ -1,7 +1,7 @@
-version = version = (1, 2, 9)
+version = version = (1, 3, 0)
 __authors__ = "Peter Burek, Johan van der Knijff, Ad de Roo"
 __version__ = '.'.join(list(map(str, version)))
-__date__ = "22/01/2024"
+__date__ = "07/03/2024"
 __copyright__ = "Copyright 2020, Lisflood Open Source"
 __maintainers__ = "Goncalo Gomes, Domenico Nappo, Valerio Lorini, Lorenzo Mentaschi, Lorenzo Alfieri"
 __status__ = "Operation"

src/lisvap/lisvapdynamic.py

@@ -135,7 +135,7 @@ def dynamic(self):
         # CM: get time for operation "Start dynamic"
         tp.timemeasure('Start dynamic')
         # CM: date corresponding to the model time step (yyyy-mm-dd hh:mm:ss)
-        self.calendar_date = self.calendar_day_start + datetime.timedelta(days=(self.currentTimeStep()) * self.DtDay)
+        self.calendar_date = self.calendar_day_start + datetime.timedelta(seconds=(self.currentTimeStep()) * self.DtSec)
         # CM: day of the year corresponding to the model time step
         self.calendar_day = int(self.calendar_date.strftime("%j"))
 
@@ -146,7 +146,7 @@ def dynamic(self):
         self.time_since_start = self.currentTimeStep() - self.firstTimeStep() + 1
 
         if settings.flags['loud']:
-            print("%-6i %10s" % (self.currentTimeStep(), self.calendar_date.strftime("%d/%m/%Y")))
+            print("%-6i %10s" % (self.currentTimeStep(), self.calendar_date.strftime("%d/%m/%Y %H:%M")))
         elif not settings.flags['checkfiles']:
             if settings.flags['quiet'] and not settings.flags['veryquiet']:
                 sys.stdout.write(".")

src/lisvap/utils/readers.py

@@ -275,9 +275,17 @@ def netcdf_step(averageyearflag, nf1, timestampflag, timestep, splitIO=False):
     # Time step, expressed as fraction of day (same as self.var.DtSec and self.var.DtDay)
     # get date of current simulation step
     current_date = calendar(timestep)
+
     if not isinstance(current_date, datetime.datetime):
-        current_date_number = date2num(begin, units=t_unit, calendar=t_cal) + ((current_date - 1) * DtDay)
-        current_date = num2date(current_date_number, t_unit, t_cal)
+        current_date_number = current_date * int(settings.binding['internal.time.frequency'])
+        init_t_unit = settings.binding['internal.time.unit']
+        init_t_cal = settings.binding['internal.time.calendar']
+        begin_date_number = date2num(begin, units=init_t_unit, calendar=init_t_cal)
+        cur_date = num2date(current_date_number, init_t_unit, init_t_cal)
+        next_date = cur_date - datetime.timedelta(seconds=DtSec)
+        cur_step = date2num(next_date, units=init_t_unit, calendar=init_t_cal)
+        current_date_number = begin_date_number + cur_step
+        current_date = num2date(current_date_number, init_t_unit, init_t_cal)
     # if reading from an average year NetCDF stack, ignore the year in current simulation date and change it to the netCDF time unit year
     if averageyearflag:
         # CM: get year from time unit in case average year is used

src/lisvap/utils/tools.py

@@ -29,6 +29,7 @@
 from pcraster.operations import scalar, defined, maptotal, ifthenelse, mapminimum, mapmaximum
 
 from . import LisfloodError
+from . import LisSettings
 from .decorators import counted
 
 
@@ -97,10 +98,16 @@ def get_calendar_configuration(netcdf_file_obj, settings=None):
     except AttributeError:  # Attribute does not exist
         t_unit = u'hours since 1990-01-01 06:00:00'
         t_cal = u'gregorian'  # Use standard calendar
+    try:
+        t_frequency = int(netcdf_file_obj.variables['time'].frequency) # get frequency
+    except AttributeError:  # Attribute does not exist
+        t_frequency = 1
     if settings is not None and not ('internal.time.unit' in settings.binding and
-                                     'internal.time.calendar' in settings.binding):
+                                     'internal.time.calendar' in settings.binding and
+                                     'internal.time.frequency' in settings.binding):
         settings.binding['internal.time.unit'] = t_unit
         settings.binding['internal.time.calendar'] = t_cal
+        settings.binding['internal.time.frequency'] = t_frequency
     return t_unit, t_cal
 
 
@@ -142,7 +149,6 @@ def date_to_int(date_in, both=False):
     the number of steps as integer is returned
     :return: number of steps as integer and input date as string
     """
-    from lisvap.utils import LisSettings
     settings = LisSettings.instance()
     binding = settings.binding
     # CM: get reference date to be used with step numbers from 'CalendarDayStart' in Settings.xml file
@@ -175,7 +181,6 @@ def checkdate(start, end):
     :param start: start date for model run (# or date as string)
     :param end: end date for model run (# or date as string)
     """
-    from . import LisSettings
     settings = LisSettings.instance()
     binding = settings.binding
     # CM: calendar date start (CalendarDayStart)

src/lisvap/utils/writers.py

@@ -63,55 +63,47 @@ def coordinates_range(start=0, nelems=1, step=1):
     return elem_array
 
 
-def get_output_parameters_monthly(start_date, timestep, current_output_index):
+def get_output_parameters_monthly(start_date, timestep, time_frequency, timestep_stride, current_output_index):
     output_index = current_output_index
-    start_yearmonth = start_date.strftime('%Y%m')
-    current_date = start_date + datetime.timedelta(days=timestep-1)
-    current_yearmonth = current_date.strftime('%Y%m')
-    filename_suffix = current_yearmonth
-    if start_yearmonth != current_yearmonth:
-        first_day_current_month = current_date.replace(day=1, hour=0, minute=0, second=0, microsecond=0)
-        current_day = current_date.replace(hour=0, minute=0, second=0, microsecond=0)
-        days_between = (current_day - first_day_current_month).days
-        if days_between == 0:
-            # Get last month because the 1st of the month belongs in the last month file
-            last_day_last_month = first_day_current_month - datetime.timedelta(days=1)
-            filename_suffix = (last_day_last_month).strftime('%Y%m')
-            num_days_last_month = last_day_last_month.day
-            # If last month was complete the idx needs to be set to the last idx of the previous file
-            if output_index > num_days_last_month:
-                output_index = num_days_last_month - 1
-        else:
-            # Since the 1st day belongs to the last year file, the remaining days need to start on index 0...
-            output_index = days_between - 1
+    current_date = start_date + datetime.timedelta(seconds=((timestep - 1) * timestep_stride))
+    filename_suffix = current_date.strftime('%Y%m')
+
+    first_date_current_month = start_date.replace(year=current_date.year, month=current_date.month)
+    seconds_between = (current_date - first_date_current_month).total_seconds()
+    num_steps_done_in_current_month = int(seconds_between / timestep_stride) + 1
+    last_date_last_month = first_date_current_month - datetime.timedelta(seconds=timestep_stride)
+    day_inside_last_month = last_date_last_month - datetime.timedelta(seconds=timestep_stride)
+
+    if current_date == first_date_current_month:
+        output_index = 0
+    elif current_date == last_date_last_month:
+        filename_suffix = day_inside_last_month.strftime('%Y%m')
+    elif current_output_index >= num_steps_done_in_current_month:
+        output_index = num_steps_done_in_current_month - 1
     return filename_suffix, output_index
 
 
-def get_output_parameters_yearly(start_date, timestep, current_output_index):
+def get_output_parameters_yearly(start_date, timestep, time_frequency, timestep_stride, current_output_index):
     output_index = current_output_index
-    start_year = start_date.strftime('%Y')
-    current_date = start_date + datetime.timedelta(days=timestep-1)
-    current_year = current_date.strftime('%Y')
-    filename_suffix = current_year
-    if start_year != current_year:
-        first_day_current_year = current_date.replace(month=1, day=1, hour=0, minute=0, second=0, microsecond=0)
-        current_day = current_date.replace(hour=0, minute=0, second=0, microsecond=0)
-        days_between = (current_day - first_day_current_year).days
-        if days_between == 0:
-            # Get last year because the 1st of January belongs in the last year file
-            last_day_last_year = first_day_current_year - datetime.timedelta(days=1)
-            filename_suffix = (last_day_last_year).strftime('%Y')
-            num_days_last_year = last_day_last_year.timetuple().tm_yday
-            # If last year was complete the idx needs to be set to the last idx of the previous file
-            if output_index > num_days_last_year:
-                output_index = num_days_last_year - 1
-        else:
-            # Since the 1st day belongs to the last year file, the remaining days need to start on index 0...
-            output_index = days_between - 1
+    current_date = start_date + datetime.timedelta(seconds=((timestep - 1) * timestep_stride))
+    filename_suffix = current_date.strftime('%Y')
+
+    first_date_current_year = start_date.replace(year=current_date.year)
+    seconds_between = (current_date - first_date_current_year).total_seconds()
+    num_steps_done_in_current_year = int(seconds_between / timestep_stride) + 1
+    last_date_last_year = first_date_current_year - datetime.timedelta(seconds=timestep_stride)
+    day_inside_last_year = last_date_last_year - datetime.timedelta(seconds=timestep_stride)
+
+    if current_date == first_date_current_year:
+        output_index = 0
+    elif current_date == last_date_last_year:
+        filename_suffix = day_inside_last_year.strftime('%Y')
+    elif current_output_index >= num_steps_done_in_current_year:
+        output_index = num_steps_done_in_current_year - 1
     return filename_suffix, output_index
 
 
-def get_output_parameters(settings, netcdf_output_file, start_date, timestep, current_output_index):
+def get_output_parameters(settings, netcdf_output_file, start_date, timestep, time_frequency, timestep_stride, current_output_index):
     output_index = current_output_index
     p = Path(netcdf_output_file)
     netfile = Path(p.parent) / Path('{}.nc'.format(p.name) if not p.name.endswith('.nc') else p.name)
@@ -120,9 +112,9 @@ def get_output_parameters(settings, netcdf_output_file, start_date, timestep, cu
     if splitIO:
         monthlyIO = settings.get_option('monthlyOutput')
         if monthlyIO:
-            filename_suffix, output_index = get_output_parameters_monthly(start_date, timestep, current_output_index)
+            filename_suffix, output_index = get_output_parameters_monthly(start_date, timestep, time_frequency, timestep_stride, current_output_index)
         else:
-            filename_suffix, output_index = get_output_parameters_yearly(start_date, timestep, current_output_index)
+            filename_suffix, output_index = get_output_parameters_yearly(start_date, timestep, time_frequency, timestep_stride, current_output_index)
         netfile = Path(p.parent) / Path('{}_{}.nc'.format(p.name, filename_suffix) if not p.name.endswith('.nc') else p.name)
     return prefix, netfile, output_index
 
@@ -140,6 +132,16 @@ def set_time_dimension(settings, netcdf_obj, time_variable_name, start_date, var
         start_date_6hourly = start_date - datetime.timedelta(hours=18)
         time.units = 'hours since %s' % start_date_6hourly.strftime('%Y-%m-%d %H:%M:%S.0')
         time.frequency = 6
+    elif 'internal.time.unit' in settings.binding:
+        internal_time_units = settings.binding['internal.time.unit']
+        # Separate the different parts of time units, ex:
+        # "seconds since 1970-01-01 00:00:00"
+        # "days since 1970-01-01 00:00:00"
+        # "hours since 1970-01-01 00:00:00"
+        # "minutes since 1970-01-01 00:00:00"
+        time_units_parts = internal_time_units.split(' ')
+        time.units = '%s since %s' % (time_units_parts[0], start_date.strftime('%Y-%m-%d %H:%M:%S.0'))
+        time.frequency = int(settings.binding['internal.time.frequency'])
     else:
         time.units = 'days since %s' % start_date.strftime('%Y-%m-%d %H:%M:%S.0')
         time.frequency = 1
@@ -276,19 +278,21 @@ def writenet(current_output_index, inputmap, netcdf_output_file, current_timeste
     netfile: output netcdf filename
     timestep:
     """
-    start_date = calendar_day_start + datetime.timedelta(days=1)
+    settings = LisSettings.instance()
+
+    timestep_stride = int(settings.binding['DtSec'])
+    time_frequency = int(settings.binding['internal.time.frequency'])
+    start_date = calendar_day_start + datetime.timedelta(seconds=timestep_stride)
     timestep = current_timestep
 
     cutmap = CutMap.instance()
     nrows = np.abs(cutmap.cuts[3] - cutmap.cuts[2])
     ncols = np.abs(cutmap.cuts[1] - cutmap.cuts[0])
 
-    settings = LisSettings.instance()
-
     time_variable = 'time'
     output6hourly = settings.get_option('output6hourly')
-    prefix, netfile, output_index = get_output_parameters(settings, netcdf_output_file, start_date,
-                                                          timestep, current_output_index)
+    prefix, netfile, output_index = get_output_parameters(settings, netcdf_output_file, start_date, timestep,
+                                                          time_frequency, timestep_stride, current_output_index)
 
     # Create and setup the netcdf file when the first map needs to be stored
     if output_index == 0 or not netfile.exists():
@@ -304,15 +308,15 @@ def writenet(current_output_index, inputmap, netcdf_output_file, current_timeste
     mapnp[np.isnan(mapnp)] = (nan_value - add_offset) * scale_factor
     if flag_time:
         # In case output6hourly==True, replicate four daily maps to get the 6 hourly output (EFCC-2316)
-        # The timestep need to increase by 4
+        # The timestep needs to increase by 4
         if output6hourly:
             time_frequency = 6
             for i in range(4):
                 map_idx = output_index * 4 + i
                 nf1.variables[time_variable][map_idx] = (timestep * 4 - 4 + i) * time_frequency
                 nf1.variables[prefix][map_idx, :, :] = mapnp
         else:  # Generate daily output
-            nf1.variables[time_variable][output_index] = timestep - 1
+            nf1.variables[time_variable][output_index] = (timestep - 1) * time_frequency  # timestep - time_frequency
             nf1.variables[prefix][output_index, :, :] = mapnp
     else:
         nf1.variables[prefix][:, :] = mapnp

src/lisvap1.py

@@ -50,10 +50,11 @@ def lisvapexe(settings):
     tp = TimeProfiler()
     step_start = settings.binding['StepStart']
     step_end = settings.binding['StepEnd']
+    timestep_stride = int(settings.binding['DtSec'])
     start_date, end_date = datetime.datetime.strptime(step_start, '%d/%m/%Y %H:%M'), datetime.datetime.strptime(step_end, '%d/%m/%Y %H:%M')
     start_date_simulation = datetime.datetime.strptime(settings.binding['CalendarDayStart'], '%d/%m/%Y %H:%M')
-    timestep_start = (start_date - start_date_simulation).days + 1
-    timestep_end = (end_date - start_date_simulation).days + 1
+    timestep_start = int((start_date - start_date_simulation).total_seconds() / timestep_stride) + 1
+    timestep_end = int((end_date - start_date_simulation).total_seconds() / timestep_stride) + 1
     checkdate('StepStart', 'StepEnd')
     print('Start date: {} ({}) - End date: {} ({})'.format(step_start, timestep_start, step_end, timestep_end))