Fix merge conflicts for flip maps. Fix maps caching unit test

src/lisflood/__init__.py

@@ -10,7 +10,7 @@
 
 __version__ = version
 __authors__ = "Ad de Roo, Emiliano Gelati, Peter Burek, Johan van der Knijff, Niko Wanders"
-__date__ = "22/06/2023"
+__date__ = "11/07/2023"
 __copyright__ = "Copyright 2019-2023, European Commission - Joint Research Centre"
 __maintainer__ = "Stefania Grimaldi, Cinzia Mazzetti, Carlo Russo, Damien Decremer, Corentin Carton De Wiart, Valerio Lorini, Ad de Roo"
 __status__ = "Operation"

src/lisflood/global_modules/add1.py

@@ -140,29 +140,28 @@ def mapattrNetCDF(name):
     filename = os.path.splitext(name)[0] + '.nc'
     nf1 = iterOpenNetcdf(filename, "Checking netcdf map \n", 'r')
     spatial_dims = ('x', 'y') if 'x' in nf1.variables else ('lon', 'lat')
-
+    x1, x2, y1, y2 = [nf1.variables[v][j] for v in spatial_dims for j in (0, 1)]
+    x_left = min(nf1.variables[spatial_dims[0]][0],nf1.variables[spatial_dims[0]][-1])
+    y_top = max(nf1.variables[spatial_dims[1]][0],nf1.variables[spatial_dims[1]][-1])
+    nf1.close()
     maskattrs = MaskAttrs.instance()
-    cell_size = maskattrs['cell']
-
-    # check cell size
-    for dim in spatial_dims:
-        if len(nf1.variables[dim]) > 1:
-            x0, x1 = [nf1.variables[dim][i] for i in (0, 1)]
-            check_cell = abs(cell_size - np.abs(x1 - x0)) # this must be same precision as pcraster.clone().cellsize()
-            if check_cell > 10**-5:
-                raise LisfloodError("Cell size different in maskmap {} and {}".format(
-                    LisSettings.instance().binding['MaskMap'], filename))
-
-    x0, y0 = [nf1.variables[dim][0] for dim in spatial_dims]
-    half_cell = cell_size / 2.
-    x = x0 - half_cell  # |
-    y = y0 + half_cell  # | coordinates of the upper left corner of the input file upper left pixel
-    cut0 = int(np.abs(maskattrs['x'] - x) / cell_size)
+
+    cell_x = np.abs(x2 - x1)
+    cell_y = np.abs(y2 - y1)
+    check_x = maskattrs['cell'] - cell_x # this must be same precision as pcraster.clone().cellsize()
+    check_y = maskattrs['cell'] - cell_y # this must be same precision as pcraster.clone().cellsize()
+
+    if abs(check_x) > 10**-5 or abs(check_y) > 10**-5:
+        raise LisfloodError("Cell size different in maskmap {} and {}".format(
+            LisSettings.instance().binding['MaskMap'], filename)
+        )
+    half_cell = maskattrs['cell'] / 2.
+    x = x_left - half_cell  # |
+    y = y_top + half_cell   # | coordinates of the upper left corner of the input file upper left pixel
+    cut0 = int(np.abs(maskattrs['x'] - x) / cell_x)
     cut1 = cut0 + maskattrs['col']
-    cut2 = int(np.abs(maskattrs['y'] - y) / cell_size)
+    cut2 = int(np.abs(maskattrs['y'] - y) / cell_y)
     cut3 = cut2 + maskattrs['row']
-
-    nf1.close()
     return cut0, cut1, cut2, cut3  # input data will be sliced using [cut0:cut1,cut2:cut3]
 
 
@@ -209,20 +208,31 @@ def loadsetclone(name):
             nf1 = iterOpenNetcdf(filename, "", "r")
             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...
+            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]
                 y1 = nf1.variables['y'][0]
+                y2 = nf1.variables['y'][-1]
             else:
                 x1 = nf1.variables['lon'][0]
                 x2 = nf1.variables['lon'][-1]
                 y1 = nf1.variables['lat'][0]
+                y2 = nf1.variables['lat'][-1]
+            x_left = min(x1,x2)
+            y_top = max(y1,y2)
 
             cell_size = np.abs(x2 - x1)/(nr_cols - 1)
-            x = x1 - cell_size / 2
-            y = y1 + cell_size / 2
+            x = x_left - cell_size / 2
+            y = y_top + cell_size / 2
             mapnp = np.array(nf1.variables[varname][0:nr_rows, 0:nr_cols])
+            # read maps using always a standard x and y reference system using x in ascending and y in descending order
+            if (y1 != y_top):   # y in in ascending order
+                warnings.warn(LisfloodWarning("Warning: map {} (binding: '{}') has y coordinates in ascending order and will be flipped vertically".format(filename, name)))
+                mapnp=np.flipud(mapnp).copy()
+            if (x1 != x_left):   # x in in descending order
+                warnings.warn(LisfloodWarning("Warning: map {} (binding: '{}') has x coordinates in descending order and will be flipped horizontally".format(filename, name)))
+                mapnp=np.fliplr(mapnp).copy()
             nf1.close()
             # setclone  row col cellsize xupleft yupleft
             setclone(nr_rows, nr_cols, cell_size, x, y)
@@ -392,9 +402,24 @@ def loadmap_base(name, pcr=False, lddflag=False, timestampflag='exact', averagey
         # Only one variable must be present in netcdf files
         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]
+
+        # read maps using always a standard x and y reference system using x in ascending and y in descending order
+        spatial_dims = ('x', 'y') if 'x' in nf1.variables else ('lon', 'lat')
+        x_flipped = nf1.variables[spatial_dims[0]][0]>nf1.variables[spatial_dims[0]][-1] 
+        y_flipped = nf1.variables[spatial_dims[1]][0]<nf1.variables[spatial_dims[1]][-1] 
+        func_y = lambda y : y
+        func_x = lambda x : x
+        # read maps using always a standard x and y reference system using x in ascending and y in descending order
+        if (y_flipped):   # y in in ascending order
+            warnings.warn(LisfloodWarning("Warning: map {} (binding: '{}') has y coordinates in ascending order and will be flipped vertically".format(filename, name)))
+            func_y = lambda y : np.flipud(y).copy()
+        if (x_flipped):   # x in in descending order
+            warnings.warn(LisfloodWarning("Warning: map {} (binding: '{}') has x coordinates in descending order and will be flipped horizontally".format(filename, name)))
+            func_x = lambda x : np.fliplr(x).copy()
+
         if not settings.timestep_init:
             # if timestep_init is missing, read netcdf as single static map
-            mapnp = nf1.variables[varname][cut2:cut3, cut0:cut1]
+            mapnp = func_x(func_y(nf1.variables[varname]))[cut2:cut3, cut0:cut1]
         else:
             if 'time' in nf1.variables:
                 # read a netcdf  (stack) - state files
@@ -456,10 +481,10 @@ def loadmap_base(name, pcr=False, lddflag=False, timestampflag='exact', averagey
                         timestepI = timestepInew
 
                 itime = np.where(nf1.variables['time'][:] == timestepI)[0][0]
-                mapnp = nf1.variables[varname][itime, cut2:cut3, cut0:cut1]
+                mapnp = func_x(func_y(nf1.variables[varname][itime]))[cut2:cut3, cut0:cut1]
             else:
                 # read a netcdf (single one)
-                mapnp = nf1.variables[varname][cut2:cut3, cut0:cut1]
+                mapnp = func_x(func_y(nf1.variables[varname]))[cut2:cut3, cut0:cut1]
 
         # masking
         try:
@@ -562,7 +587,21 @@ def loadLAI(value, pcrvalue, i, pcr=False):
         # Only one variable must be present in netcdf files
         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]
-        mapnp = nf1.variables[varname][i, cut2:cut3, cut0:cut1]
+
+        # read maps using always a standard x and y reference system using x in ascending and y in descending order
+        spatial_dims = ('x', 'y') if 'x' in nf1.variables else ('lon', 'lat')
+        x_flipped = nf1.variables[spatial_dims[0]][0]>nf1.variables[spatial_dims[0]][-1] 
+        y_flipped = nf1.variables[spatial_dims[1]][0]<nf1.variables[spatial_dims[1]][-1] 
+        func_y = lambda y : y
+        func_x = lambda x : x
+        # read maps using always a standard x and y reference system using x in ascending and y in descending order
+        if (y_flipped):   # y in in ascending order
+            warnings.warn(LisfloodWarning("Warning: map {} ({} maps) has y coordinates in ascending order and will be flipped vertically".format(filename, 'LAI')))
+            func_y = lambda y : np.flipud(y).copy()
+        if (x_flipped):   # x in in descending order
+            warnings.warn(LisfloodWarning("Warning: map {} ({} maps) has x coordinates in descending order and will be flipped horizontally".format(filename, 'LAI')))
+            func_x = lambda x : np.fliplr(x).copy()
+        mapnp = func_x(func_y(nf1.variables[varname][i]))[cut2:cut3, cut0:cut1]
         nf1.close()
         mapC = compressArray(mapnp, pcr=False, name=filename)
         # mapnp[np.isnan(mapnp)] = -9999
@@ -695,9 +734,23 @@ def readnetcdf(name, time, timestampflag='exact', averageyearflag=False):
     # get index of timestep in netCDF file corresponding to current simulation date
     current_ncdf_index = np.where(t_steps == current_ncdf_step)[0][0]
 
+    # read maps using always a standard x and y reference system using x in ascending and y in descending order
+    spatial_dims = ('x', 'y') if 'x' in nf1.variables else ('lon', 'lat')
+    x_flipped = nf1.variables[spatial_dims[0]][0]>nf1.variables[spatial_dims[0]][-1] 
+    y_flipped = nf1.variables[spatial_dims[1]][0]<nf1.variables[spatial_dims[1]][-1] 
+    func_y = lambda y : y
+    func_x = lambda x : x
+    # read maps using always a standard x and y reference system using x in ascending and y in descending order
+    if (y_flipped):   # y in in ascending order
+        warnings.warn(LisfloodWarning("Warning: map {} (binding: '{}') has y coordinates in ascending order and will be flipped vertically".format(filename, name)))
+        func_y = lambda y : np.flipud(y).copy()
+    if (x_flipped):   # x in in descending order
+        warnings.warn(LisfloodWarning("Warning: map {} (binding: '{}') has x coordinates in descending order and will be flipped horizontally".format(filename, name)))
+        func_x = lambda x : np.fliplr(x).copy()
+
     # crop map if subcatchment
     cutmap = CutMap.instance()
-    mapnp = nf1.variables[variable_name][current_ncdf_index, cutmap.cuts[2]:cutmap.cuts[3], cutmap.cuts[0]:cutmap.cuts[1]]
+    mapnp = func_x(func_y(nf1.variables[variable_name][current_ncdf_index]))[cutmap.cuts[2]:cutmap.cuts[3], cutmap.cuts[0]:cutmap.cuts[1]]
     nf1.close()
 
     mapC = compressArray(mapnp, pcr=False, name=filename)

src/lisflood/global_modules/netcdf.py

@@ -21,8 +21,8 @@
 
 from .. import __authors__, __version__, __date__, __status__, __institution__
 
-def mask_array_np(data, mask, crop, name, valid_min, valid_max):
-    data_cut = data[:, crop[2]:crop[3], crop[0]:crop[1]]
+def mask_array_np(data, mask, crop, name, valid_min, valid_max, func_x, func_y):
+    data_cut = func_x(func_y(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)))
@@ -41,7 +41,7 @@ def mask_array_np(data, mask, crop, name, valid_min, valid_max):
 
 
 
-def mask_array(data, mask, crop, core_dims, name, valid_min, valid_max):
+def mask_array(data, mask, crop, core_dims, name, valid_min, valid_max, func_x, func_y):
     n_data = int(mask.sum())
     masked_data = xr.apply_ufunc(mask_array_np, data,
                                  dask='parallelized',
@@ -50,13 +50,17 @@ def mask_array(data, mask, crop, core_dims, name, valid_min, valid_max):
                                  output_dtypes=[data.dtype],
                                  output_core_dims=[['z']],
                                  dask_gufunc_kwargs = dict(output_sizes={'z': n_data}),
-                                 kwargs={'mask': mask, 'crop': crop, 'name': name, 'valid_min': valid_min, 'valid_max': valid_max})
+                                 kwargs={'mask': mask, 'crop': crop, 'name': name, 
+                                         'valid_min': valid_min, 'valid_max': valid_max,
+                                         'func_x': func_x, 'func_y':func_y})
     return masked_data
 
 
-def compress_xarray(mask, crop, data, name, valid_min, valid_max):
+def compress_xarray(mask, crop, data, name, valid_min, valid_max, func_x, func_y):
     core_dims = get_core_dims(data.dims)
-    masked_data = mask_array(data, mask, crop, core_dims=core_dims, name=name, valid_min=valid_min, valid_max=valid_max)
+    masked_data = mask_array(data, mask, crop, core_dims=core_dims, name=name, 
+                             valid_min=valid_min, valid_max=valid_max,
+                             func_x=func_x, func_y=func_y)
     return masked_data
 
 
@@ -203,6 +207,20 @@ def __init__(self, data_path, time_chunk, dates, indexer=None, climatology=False
 
         self.index_map = map_dates_index(dates, da.time, indexer, climatology)
 
+        # read maps using always a standard x and y reference system using x in ascending and y in descending order
+        spatial_dims = ('x', 'y') if 'x' in ds.variables else ('lon', 'lat')
+        x_flipped = ds.variables[spatial_dims[0]][0]>ds.variables[spatial_dims[0]][-1] 
+        y_flipped = ds.variables[spatial_dims[1]][0]<ds.variables[spatial_dims[1]][-1] 
+        func_y = lambda y : y
+        func_x = lambda x : x
+        # read maps using always a standard x and y reference system using x in ascending and y in descending order
+        if (y_flipped):   # y in in ascending order
+            warnings.warn(LisfloodWarning("Warning: map {} (var_name: '{}') has y coordinates in ascending order and will be flipped vertically".format(data_path, var_name)))
+            func_y = lambda y : np.flipud(y).copy()
+        if (x_flipped):   # x in in descending order
+            warnings.warn(LisfloodWarning("Warning: map {} (var_name: '{}') has x coordinates in descending order and will be flipped horizontally".format(data_path, var_name)))
+            func_x = lambda x : np.fliplr(x).copy()
+
         # compress dataset (remove missing values and flatten the array)
         maskinfo = MaskInfo.instance()
         mask = np.logical_not(maskinfo.info.mask)
@@ -219,7 +237,7 @@ def __init__(self, data_path, time_chunk, dates, indexer=None, climatology=False
         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
+        self.dataset = compress_xarray(mask, crop, da, var_name, valid_min_scaled, valid_max_scaled, func_x, func_y) # final dataset to store
 
         # initialise class variables and load first chunk
         self.init_chunks(self.dataset, time_chunk)
@@ -343,6 +361,19 @@ def read_lat_from_template_cached(netcdf_template, proj4_params):
 def read_lat_from_template_base(netcdf_template, proj4_params):
     nc_template = netcdf_template + ".nc" if not netcdf_template.endswith('.nc') else netcdf_template
     with xr.open_dataset(nc_template) as nc:
+        # read maps using always a standard x and y reference system using x in ascending and y in descending order
+        spatial_dims = ('x', 'y') if 'x' in nc.dims else ('lon', 'lat')
+        x_flipped = nc.variables[spatial_dims[0]][0]>nc.variables[spatial_dims[0]][-1] 
+        y_flipped = nc.variables[spatial_dims[1]][0]<nc.variables[spatial_dims[1]][-1] 
+        func_y = lambda y : y
+        func_x = lambda x : x
+        # read maps using always a standard x and y reference system using x in ascending and y in descending order
+        if (y_flipped):   # y in in ascending order
+            warnings.warn(LisfloodWarning("Warning: map {} (binding: '{}') has y coordinates in ascending order and will be flipped vertically".format(nc_template,'netCDFtemplate')))
+            func_y = lambda y : np.flip(y).copy()
+        if (x_flipped):   # x in in descending order
+            warnings.warn(LisfloodWarning("Warning: map {} (binding: '{}') has x coordinates in descending order and will be flipped horizontally".format(nc_template, 'netCDFtemplate')))
+            func_x = lambda x : np.flip(x).copy()
         if all([co in nc.dims for co in ("x", "y")]):
             try:
                 # look for the projection variable
@@ -357,9 +388,10 @@ def read_lat_from_template_base(netcdf_template, proj4_params):
 
             # projection object obtained from the PROJ4 string
             projection = Proj(proj4_params)
-            _, lat_deg = projection(*coordinatesLand(nc.x.values, nc.y.values), inverse=True)  # latitude (degrees)
+
+            _, lat_deg = projection(*coordinatesLand(func_x(nc.x.values), func_y(nc.y.values)), inverse=True)  # latitude (degrees)
         else:
-            _, lat_deg = coordinatesLand(nc.lon.values, nc.lat.values)  # latitude (degrees)
+            _, lat_deg = coordinatesLand(func_x(nc.lon.values), func_y(nc.lat.values))  # latitude (degrees)
 
     return lat_deg
 
@@ -475,7 +507,8 @@ def write_netcdf_header(settings,
     nrow = np.abs(cutmap.cuts[3] - cutmap.cuts[2])
     ncol = np.abs(cutmap.cuts[1] - cutmap.cuts[0])
 
-    dim_lat_y, dim_lon_x = get_core_dims(meta_netcdf.data)
+    # since meta_netcdf is the in correct x and y reference system, variables will be correctly cut and written in the same order (x ascending and y descending)
+    dim_lat_y, dim_lon_x = get_core_dims(meta_netcdf.data)    
     latlon_coords = get_space_coords(meta_netcdf, dim_lat_y, dim_lon_x)
 
     if dim_lon_x in meta_netcdf.data:

src/lisflood/global_modules/settings.py

@@ -297,8 +297,19 @@ def __init__(self, uid):
                 self.data[var] = {k: v for k, v in iteritems(nf1.variables[var].__dict__) if k != '_FillValue'}
             core_dims = get_core_dims(self.data)
             self.coords = {}
-            for dim in core_dims:
-                self.coords[dim] = nf1.variables[dim][:]
+            x_flipped = nf1.variables[core_dims[1]][0]>nf1.variables[core_dims[1]][-1] 
+            y_flipped = nf1.variables[core_dims[0]][0]<nf1.variables[core_dims[0]][-1] 
+            func_y = lambda y : y
+            func_x = lambda x : x
+            # write maps using always a standard x and y reference system using x in ascending and y in descending order
+            if (y_flipped):   # y in in ascending order
+                warnings.warn(LisfloodWarning("Warning: map {} (netCDFtemplate) has y coordinates in ascending order and will be flipped vertically".format(filename)))
+                func_y = lambda y : np.flipud(y).copy()
+            if (x_flipped):   # x in in descending order
+                warnings.warn(LisfloodWarning("Warning: map {} (netCDFtemplate) has x coordinates in descending order and will be flipped horizontally".format(filename)))
+                func_x = lambda x : np.fliplr(x).copy()
+            self.coords[core_dims[1]] = func_x(nf1.variables[core_dims[1]][:])
+            self.coords[core_dims[0]] = func_y(nf1.variables[core_dims[0]][:])
             nf1.close()
             return
         except (KeyError, IOError, IndexError, Exception):

src/lisfloodSettings_reference.xml

@@ -1725,6 +1725,7 @@ This is necessary when using projected coordinates (x,y)
 </textvar>
 
 <textvar name="OutputMapsChunks" value="$(OutputMapsChunks)"/>
+<textvar name="OutputMapsDataType" value="$(OutputMapsDataType)"/>
 <textvar name="NetCDFTimeChunks" value="$(NetCDFTimeChunks)"/>
 <textvar name="MapsCaching" value="$(MapsCaching)"/>