diff --git a/VERSION b/VERSION index cd9b8f55..8c7fafd3 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -4.1.2 \ No newline at end of file +4.1.3 \ No newline at end of file diff --git a/requirements.txt b/requirements.txt index 1a350afa..12d00c95 100644 --- a/requirements.txt +++ b/requirements.txt @@ -18,7 +18,7 @@ numba>=0.54.0 numexpr>=2.7.0 numpy>=1.21.0 packaging>=19.2 -pandas>=0.25.1 +pandas>=1 pathlib2>=2.3.5 pkginfo>=1.5.0.1 pluggy>=0.13.0 diff --git a/src/lisflood/Lisflood_dynamic.py b/src/lisflood/Lisflood_dynamic.py index 364f0a40..424ac81c 100644 --- a/src/lisflood/Lisflood_dynamic.py +++ b/src/lisflood/Lisflood_dynamic.py @@ -247,7 +247,7 @@ def splitlanduse(array1, array2=None, array3=None): ftemp1 = open(nomefile, 'w+') nelements = len(self.ChanM3) for i in range(0,nelements-1): - if hasattr(self,'CrossSection2Area'): + if hasattr(self,'CrossSection2Area') and not(option['InitLisflood']): print(i, self.TotalCrossSectionArea[i], self.CrossSection2Area[i], self.ChanM3[i], \ self.Chan2M3Kin[i], file=ftemp1) else: diff --git a/src/lisflood/Lisflood_initial.py b/src/lisflood/Lisflood_initial.py index 0a77589f..a199a6cd 100644 --- a/src/lisflood/Lisflood_initial.py +++ b/src/lisflood/Lisflood_initial.py @@ -238,7 +238,7 @@ def __init__(self): ftemp1 = open(nomefile, 'w+') nelements = len(self.ChanM3) for i in range(0, nelements - 1): - if hasattr(self, 'CrossSection2Area'): + if hasattr(self, 'CrossSection2Area') and not(option['InitLisflood']): print(i, self.TotalCrossSectionArea[i], self.CrossSection2Area[i], self.ChanM3[i], self.Chan2M3Kin[i], file=ftemp1) else: print(i, self.TotalCrossSectionArea[i], self.ChanM3[i], file=ftemp1) diff --git a/src/lisflood/global_modules/add1.py b/src/lisflood/global_modules/add1.py index bbd3efd2..11162ea3 100755 --- a/src/lisflood/global_modules/add1.py +++ b/src/lisflood/global_modules/add1.py @@ -204,8 +204,9 @@ def loadsetclone(name): # try to read a netcdf file filename = os.path.splitext(binding[name])[0] + '.nc' nf1 = iterOpenNetcdf(filename, "", "r") - value = listitems(nf1.variables)[-1][0] # get the last variable name - nr_rows, nr_cols = nf1.variables[value].shape # just use shape to know rows and cols... + num_dims = 3 if 'time' in nf1.variables else 2 + varname = [v for v in nf1.variables if len(nf1.variables[v].dimensions) == num_dims][0] + nr_rows, nr_cols = nf1.variables[varname].shape # just use shape to know rows and cols... if 'x' in nf1.variables: x1 = nf1.variables['x'][0] x2 = nf1.variables['x'][-1] @@ -218,7 +219,7 @@ def loadsetclone(name): cell_size = np.abs(x2 - x1)/(nr_cols - 1) x = x1 - cell_size / 2 y = y1 + cell_size / 2 - mapnp = np.array(nf1.variables[value][0:nr_rows, 0:nr_cols]) + mapnp = np.array(nf1.variables[varname][0:nr_rows, 0:nr_cols]) nf1.close() # setclone row col cellsize xupleft yupleft setclone(nr_rows, nr_cols, cell_size, x, y) diff --git a/src/lisflood/global_modules/default_options.py b/src/lisflood/global_modules/default_options.py index 7126c023..7ae6c4c7 100644 --- a/src/lisflood/global_modules/default_options.py +++ b/src/lisflood/global_modules/default_options.py @@ -897,126 +897,126 @@ 'TavgMapsOut': ReportedMap(name='TavgMapsOut', output_var='Tavg', unit='degree', end=[], steps=[], all=['repTavgMaps'], restrictoption=[], monthly=False, yearly=False), - 'Theta1End': ReportedMap(name='Theta1End', output_var='Theta1a[0]', unit='mm', + 'Theta1End': ReportedMap(name='Theta1End', output_var='Theta1a[0]', unit='', end=['repEndMaps'], steps=[], all=[], restrictoption=[], monthly=False, yearly=False), 'Theta1ForestEnd': ReportedMap(name='Theta1ForestEnd', output_var='Theta1a[1]', - unit='mm', end=['repEndMaps'], steps=[], + unit='', end=['repEndMaps'], steps=[], all=[], restrictoption=[], monthly=False, yearly=False), 'Theta1ForestMaps': ReportedMap(name='Theta1ForestMaps', output_var='Theta1a[1]', - unit='mm', end=[], steps=[], + unit='', end=[], steps=[], all=['repThetaForestMaps','repThetaMaps'], restrictoption=['nonInit'], monthly=False, yearly=False), 'Theta1ForestState': ReportedMap(name='Theta1ForestState', - output_var='Theta1a[1]', unit='mm', end=[], + output_var='Theta1a[1]', unit='', end=[], steps=['repStateMaps'], all=[], restrictoption=['nonInit'], monthly=False, yearly=False), 'Theta1IrrigationEnd': ReportedMap(name='Theta1IrrigationEnd', - output_var='Theta1a[2]', unit='mm', + output_var='Theta1a[2]', unit='', end=['repEndMaps'], steps=[], all=[], restrictoption=[], monthly=False, yearly=False), 'Theta1IrrigationMaps': ReportedMap(name='Theta1IrrigationMaps', - output_var='Theta1a[2]', unit='mm', + output_var='Theta1a[2]', unit='', end=[], steps=[], all=['repThetaIrrigationMaps','repThetaMaps'], restrictoption=['nonInit'], monthly=False, yearly=False), 'Theta1IrrigationState': ReportedMap(name='Theta1IrrigationState', - output_var='Theta1a[2]', unit='mm', + output_var='Theta1a[2]', unit='', end=[], steps=['repStateMaps'], all=[], restrictoption=['nonInit'], monthly=False, yearly=False), - 'Theta1Maps': ReportedMap(name='Theta1Maps', output_var='Theta1a[0]', unit='mm', + 'Theta1Maps': ReportedMap(name='Theta1Maps', output_var='Theta1a[0]', unit='', end=[], steps=[], all=['repThetaMaps','repE2O2'], restrictoption=['nonInit'], monthly=False, yearly=False), 'Theta1State': ReportedMap(name='Theta1State', output_var='Theta1a[0]', - unit='mm', end=[], steps=['repStateMaps'], + unit='', end=[], steps=['repStateMaps'], all=[], restrictoption=['nonInit'], monthly=False, yearly=False), - 'Theta2End': ReportedMap(name='Theta2End', output_var='Theta1b[0]', unit='mm', + 'Theta2End': ReportedMap(name='Theta2End', output_var='Theta1b[0]', unit='', end=['repEndMaps'], steps=[], all=[], restrictoption=[], monthly=False, yearly=False), 'Theta2ForestEnd': ReportedMap(name='Theta2ForestEnd', output_var='Theta1b[1]', - unit='mm', end=['repEndMaps'], steps=[], + unit='', end=['repEndMaps'], steps=[], all=[], restrictoption=[], monthly=False, yearly=False), 'Theta2ForestMaps': ReportedMap(name='Theta2ForestMaps', output_var='Theta1b[1]', - unit='mm', end=[], steps=[], + unit='', end=[], steps=[], all=['repThetaForestMaps','repThetaMaps'], restrictoption=['nonInit'], monthly=False, yearly=False), 'Theta2ForestState': ReportedMap(name='Theta2ForestState', - output_var='Theta1b[1]', unit='mm', end=[], + output_var='Theta1b[1]', unit='', end=[], steps=['repStateMaps'], all=[], restrictoption=['nonInit'], monthly=False, yearly=False), 'Theta2IrrigationEnd': ReportedMap(name='Theta2IrrigationEnd', - output_var='Theta1b[2]', unit='mm', + output_var='Theta1b[2]', unit='', end=['repEndMaps'], steps=[], all=[], restrictoption=[], monthly=False, yearly=False), 'Theta2IrrigationMaps': ReportedMap(name='Theta2IrrigationMaps', - output_var='Theta1b[2]', unit='mm', + output_var='Theta1b[2]', unit='', end=[], steps=[], all=['repThetaIrrigationMaps','repThetaMaps'], restrictoption=['nonInit'], monthly=False, yearly=False), 'Theta2IrrigationState': ReportedMap(name='Theta2IrrigationState', - output_var='Theta1b[2]', unit='mm', + output_var='Theta1b[2]', unit='', end=[], steps=['repStateMaps'], all=[], restrictoption=['nonInit'], monthly=False, yearly=False), - 'Theta2Maps': ReportedMap(name='Theta2Maps', output_var='Theta1b[0]', unit='mm', + 'Theta2Maps': ReportedMap(name='Theta2Maps', output_var='Theta1b[0]', unit='', end=[], steps=[], all=['repThetaMaps','repE2O2'], restrictoption=['nonInit'], monthly=False, yearly=False), 'Theta2State': ReportedMap(name='Theta2State', output_var='Theta1b[0]', - unit='mm', end=[], steps=['repStateMaps'], + unit='', end=[], steps=['repStateMaps'], all=[], restrictoption=['nonInit'], monthly=False, yearly=False), - 'Theta3End': ReportedMap(name='Theta3End', output_var='Theta2[0]', unit='mm', + 'Theta3End': ReportedMap(name='Theta3End', output_var='Theta2[0]', unit='', end=['repEndMaps'], steps=[], all=[], restrictoption=[], monthly=False, yearly=False), 'Theta3ForestEnd': ReportedMap(name='Theta3ForestEnd', output_var='Theta2[1]', - unit='mm', end=['repEndMaps'], steps=[], + unit='', end=['repEndMaps'], steps=[], all=[], restrictoption=[], monthly=False, yearly=False), 'Theta3ForestMaps': ReportedMap(name='Theta3ForestMaps', output_var='Theta2[1]', - unit='mm', end=[], steps=[], + unit='', end=[], steps=[], all=['repThetaForestMaps','repThetaMaps'], restrictoption=['nonInit'], monthly=False, yearly=False), 'Theta3ForestState': ReportedMap(name='Theta3ForestState', - output_var='Theta2[1]', unit='mm', end=[], + output_var='Theta2[1]', unit='', end=[], steps=['repStateMaps'], all=[], restrictoption=['nonInit'], monthly=False, yearly=False), 'Theta3IrrigationEnd': ReportedMap(name='Theta3IrrigationEnd', - output_var='Theta2[2]', unit='mm', + output_var='Theta2[2]', unit='', end=['repEndMaps'], steps=[], all=[], restrictoption=[], monthly=False, yearly=False), 'Theta3IrrigationMaps': ReportedMap(name='Theta3IrrigationMaps', - output_var='Theta2[2]', unit='mm', end=[], + output_var='Theta2[2]', unit='', end=[], steps=[], all=['repThetaIrrigationMaps','repThetaMaps'], restrictoption=['nonInit'], monthly=False, yearly=False), 'Theta3IrrigationState': ReportedMap(name='Theta3IrrigationState', - output_var='Theta2[2]', unit='mm', + output_var='Theta2[2]', unit='', end=[], steps=['repStateMaps'], all=[], restrictoption=['nonInit'], monthly=False, yearly=False), - 'Theta3Maps': ReportedMap(name='Theta3Maps', output_var='Theta2[0]', unit='mm', + 'Theta3Maps': ReportedMap(name='Theta3Maps', output_var='Theta2[0]', unit='', end=[], steps=[], all=['repThetaMaps','repE2O2'], restrictoption=['nonInit'], monthly=False, yearly=False), 'Theta3State': ReportedMap(name='Theta3State', output_var='Theta2[0]', - unit='mm', end=[], steps=['repStateMaps'], + unit='', end=[], steps=['repStateMaps'], all=[], restrictoption=['nonInit'], monthly=False, yearly=False), 'TotalAbsGroundwater': ReportedMap(name='TotalAbsGroundwater', diff --git a/src/lisflood/global_modules/netcdf.py b/src/lisflood/global_modules/netcdf.py index 9a53ebb9..a453f6b1 100644 --- a/src/lisflood/global_modules/netcdf.py +++ b/src/lisflood/global_modules/netcdf.py @@ -1,5 +1,6 @@ import os import glob +import warnings import xarray as xr import numpy as np import datetime @@ -16,15 +17,30 @@ from .zusatz import iterOpenNetcdf # from .add1 import * from .add1 import nanCheckMap, decompress +from netCDF4 import default_fillvals -def mask_array_np(data, mask, crop): +def mask_array_np(data, mask, crop, name, valid_min, valid_max): data_cut = data[:, crop[2]:crop[3], crop[0]:crop[1]] + if (valid_min is not None): + if (data_cut < valid_min).sum() > 0: + warnings.warn(LisfloodWarning('Data in var "{}" contains values out of valid range (valid_min)'.format(name))) + if np.issubdtype(data.dtype, np.floating): + data_cut[(data_cut < valid_min)] = np.nan + else: + data_cut[(data_cut < valid_min)] = default_fillvals[data_cut.dtype.str[1:]] + if (valid_max is not None): + if (data_cut > valid_max).sum() > 0: + warnings.warn(LisfloodWarning('Data in var "{}" contains values out of valid range (valid_max)'.format(name))) + if np.issubdtype(data.dtype, np.floating): + data_cut[(data_cut > valid_max)] = np.nan + else: + data_cut[(data_cut > valid_max)] = default_fillvals[data_cut.dtype.str[1:]] return data_cut[:, mask] -def mask_array(data, mask, crop, core_dims): +def mask_array(data, mask, crop, core_dims, name, valid_min, valid_max): n_data = int(mask.sum()) masked_data = xr.apply_ufunc(mask_array_np, data, dask='parallelized', @@ -33,13 +49,13 @@ def mask_array(data, mask, crop, core_dims): output_dtypes=[data.dtype], output_core_dims=[['z']], dask_gufunc_kwargs = dict(output_sizes={'z': n_data}), - kwargs={'mask': mask, 'crop': crop}) + kwargs={'mask': mask, 'crop': crop, 'name': name, 'valid_min': valid_min, 'valid_max': valid_max}) return masked_data -def compress_xarray(mask, crop, data): +def compress_xarray(mask, crop, data, name, valid_min, valid_max): core_dims = get_core_dims(data.dims) - masked_data = mask_array(data, mask, crop, core_dims=core_dims) + masked_data = mask_array(data, mask, crop, core_dims=core_dims, name=name, valid_min=valid_min, valid_max=valid_max) return masked_data @@ -154,7 +170,8 @@ def __init__(self, data_path, time_chunk, dates, indexer=None, climatology=False if time_chunk != 'auto' and time_chunk is not None: time_chunk = int(time_chunk) data_path = data_path + ".nc" if not data_path.endswith('.nc') else data_path - ds = xr.open_mfdataset(data_path, engine='netcdf4', chunks={'time': time_chunk}, combine='by_coords') + ds = xr.open_mfdataset(data_path, engine='netcdf4', chunks={'time': time_chunk}, combine='by_coords', + mask_and_scale=True) # check calendar type check_dataset_calendar_type(ds, data_path) @@ -174,7 +191,18 @@ def __init__(self, data_path, time_chunk, dates, indexer=None, climatology=False mask = np.logical_not(maskinfo.info.mask) cutmap = CutMap.instance() crop = cutmap.cuts - self.dataset = compress_xarray(mask, crop, da) # final dataset to store + + # in case the dataset contains valid_min and valid_max set, store here the scaled adn offset values of them + scale_factor=da.encoding['scale_factor'] if 'scale_factor' in da.encoding else 1 + add_offset=da.encoding['add_offset'] if 'add_offset' in da.encoding else 0 + valid_min_scaled = None + valid_max_scaled = None + settings = LisSettings.instance() + flags = settings.flags + if not flags['skipvalreplace']: + valid_min_scaled = (da.attrs['valid_min']*scale_factor+add_offset) if 'valid_min' in da.attrs else None + valid_max_scaled = (da.attrs['valid_max']*scale_factor+add_offset) if 'valid_max' in da.attrs else None + self.dataset = compress_xarray(mask, crop, da, var_name, valid_min_scaled, valid_max_scaled) # final dataset to store # initialise class variables and load first chunk self.init_chunks(self.dataset, time_chunk) @@ -213,9 +241,18 @@ def __getitem__(self, step): if local_index < 0: msg = 'local step cannot be negative! step: {}, chunk: {} - {}', local_index, self.chunk_index[self.ichunk], self.chunk_index[self.ichunk+1] - LisfloodError(msg) + raise LisfloodError(msg) data = self.dataset_chunk.values[local_index] + if np.issubdtype(data.dtype, np.floating): + if (np.isnan(data).any()): + #warnings.warn(LisfloodWarning('Data in var "{}" contains NaN values or values out of valid range inside mask map for step: {}'.format(self.dataset.name,local_index))) + raise LisfloodError('Data in var "{}" contains NaN values or values out of valid range inside mask map for step: {}'.format(self.dataset.name,local_index)) + else: + if (data==default_fillvals[data.dtype.str[1:]]).any(): + #warnings.warn(LisfloodWarning('Data in var "{}" contains NaN values or values out of valid range inside mask map for step: {}'.format(self.dataset.name,local_index))) + raise LisfloodError('Data in var "{}" contains missing values or values out of valid range inside mask map for step: {}'.format(self.dataset.name,local_index)) + return data diff --git a/src/lisflood/global_modules/settings.py b/src/lisflood/global_modules/settings.py index 994d2aee..e374463d 100755 --- a/src/lisflood/global_modules/settings.py +++ b/src/lisflood/global_modules/settings.py @@ -36,7 +36,7 @@ import xml.dom.minidom import pcraster from netCDF4 import Dataset, date2num, num2date -from pandas.core.tools.datetimes import parsing +from pandas import to_datetime import numpy as np from .errors import LisfloodError, LisfloodWarning, LisfloodFileError @@ -197,7 +197,7 @@ def get_core_dims(dims): core_dims = ('lat', 'lon') else: msg = 'Core dimension in netcdf file not recognised! Expecting (y, x) or (lat, lon), have '+str(dims) - LisfloodError(msg) + raise LisfloodError(msg) return core_dims @@ -331,7 +331,7 @@ def _check_timestep_init(self): float(self.timestep_init) except ValueError: try: - parsing.parse_time_string(self.timestep_init, dayfirst=True) + to_datetime(self.timestep_init, dayfirst=True).to_pydatetime() except ValueError: raise LisfloodError('Option timestepInit was not parsable. Must be integer or date string: {}'.format(self.timestep_init)) else: @@ -411,12 +411,13 @@ def _out_dir(user_settings): def _flags(sys_args): flags = OrderedDict([('quiet', False), ('veryquiet', False), ('loud', False), ('checkfiles', False), ('noheader', False), ('printtime', False), - ('debug', False), ('nancheck', False), ('initonly', False)]) + ('debug', False), ('nancheck', False), ('initonly', False), + ('skipvalreplace', False)]) @cached def _flags(argz): try: - opts, arguments = getopt.getopt(argz, 'qvlchtdni', list(flags.keys())) + opts, arguments = getopt.getopt(argz, 'qvlchtdnis', list(flags.keys())) except getopt.GetoptError: from ..main import usage usage() @@ -426,7 +427,7 @@ def _flags(argz): ('-l', '--loud'), ('-c', '--checkfiles'), ('-h', '--noheader'), ('-t', '--printtime'), ('-d', '--debug'), ('-n', '--nancheck'), - ('-i', '--initonly')): + ('-i', '--initonly'), ('-s', '--skipvalreplace')): if o in opt: flags[opt[1].lstrip('--')] = True break @@ -616,7 +617,7 @@ def calendar(date_in, calendar_type='proleptic_gregorian'): # try reading a date in one of available formats try: _t_units = "hours since 1970-01-01 00:00:00" # units used for date type conversion (datetime.datetime -> calendar-specific if needed) - date = parsing.parse_time_string(date_in, dayfirst=True)[0] # datetime.datetime type + date = to_datetime(date_in, dayfirst=True).to_pydatetime() # datetime.datetime type step = date2num(date, _t_units, calendar_type) # float type return num2date(step, _t_units, calendar_type) # calendar-dependent type from netCDF4.netcdftime._netcdftime module except: diff --git a/src/lisflood/main.py b/src/lisflood/main.py index c95954de..bfff7425 100755 --- a/src/lisflood/main.py +++ b/src/lisflood/main.py @@ -173,13 +173,15 @@ def usage(): Arguments list: settings.xml settings file - -q --quiet output progression given as . - -v --veryquiet no output progression is given - -l --loud output progression given as time step, date and discharge - -c --check input maps and stack maps are checked, output for each input map BUT no model run - -h --noheader .tss file have no header and start immediately with the time series - -d --debug debug outputs - -i --initonly only run initialisation, not the dynamic loop + -q --quiet output progression given as . + -v --veryquiet no output progression is given + -l --loud output progression given as time step, date and discharge + -c --checkfiles input maps and stack maps are checked, output for each input map BUT no model run + -n --nancheck check NaN values in output maps + -h --noheader .tss file have no header and start immediately with the time series + -d --debug debug outputs + -i --initonly only run initialisation, not the dynamic loop + -s --skipvalreplace skip replacement of invalid values in meteo input maps (ignore valid_min and valid_max) """) sys.exit(1)