improved definition coordinate grids

Author: StefaniaGrimaldi

Diff for: src/lisfloodutilities/readers/nc.py

@@ -19,7 +19,7 @@ def mv(self):
         for variable in self.ds.variables.values():
             if len(variable.dims) < 2:
                 continue
-            return variable.attrs['_FillValue'] if '_FillValue' in variable.attrs else variable.attrs['missing_value']
+            return variable.attrs['_FillValue'] if '_FillValue' in variable.attrs else -9999
 
     @property
     def data(self):

Diff for: src/lisfloodutilities/waterregions/define_waterregions.py

@@ -129,7 +129,7 @@ def define_waterregions(calib_points=None, countries_id=None, ldd=None, waterreg
           os.remove(ldd_pcr)
        except:
           pass
-       convnc2pcr(ldd,ldd_pcr,is_ldd=True) 
+       convnc2pcr(ldd,ldd_pcr,is_ldd=True)  
     else:
        ldd_pcr=ldd   
     if countries_id[-3:]=='.nc':

Diff for: src/lisfloodutilities/waterregions/verify_waterregions.py

@@ -52,19 +52,44 @@ def verify_waterregions(calib_catchments=None, waterregions=None):
        waterregions_map = Dataset(waterregions,'r','format=NETCDF4_classic')
 
        for v in waterregions_map.variables:
+           if (v == 'y' or v == 'lat'):
+              ywr = waterregions_map.variables[v][:]
+              wr_yresolution = np.abs((ywr[-1] - ywr[0])/(len(ywr[:])-1.0))
+           if (v == 'x' or v == 'lon'):
+              xwr = waterregions_map.variables[v][:]
+              wr_xresolution = np.abs((xwr[-1] - xwr[0])/(len(xwr[:])-1.0))
            if len(waterregions_map.variables[v].dimensions) == 2:
-              wr = waterregions_map.variables[v][:]             
+              wr = waterregions_map.variables[v][:]          
 
+       flipud = 0
        calibration_catchments_map = Dataset(calib_catchments,'r','format=NETCDF4_classic')
        for v in calibration_catchments_map.variables:
+           if (v == 'y' or v == 'lat'):
+              ycc = calibration_catchments_map.variables[v][:]
+              cc_resolution =np.abs((ycc[-1] - ywr[0])/(len(ycc[:])-1.0))
+              if round(ycc[0],1) != round(ywr[0],1):
+                 flipud = 1
+                 check_grid_y = np.amax(np.abs(np.flip(ycc)-ywr))
+                 if (check_grid_y > wr_yresolution):
+                      print('Error: the two maps do not have the same coordinates. Please, check the input maps.')
+                      exit()
+           if (v == 'x' or v == 'lon'):
+              xcc = calibration_catchments_map.variables[v][:]
+              check_grid_x = np.amax(np.abs(xcc-xwr))
+              if (check_grid_x > wr_xresolution):
+                      print('Error: the two maps do not have the same coordinates. Please, check the input maps.')
+                      exit()                     
            if len(calibration_catchments_map.variables[v].dimensions) == 2:
               cal_catch = calibration_catchments_map.variables[v][:]
+              if flipud == 1:
+                 print('Warning: one of the maps has inverted y-axis, is this an intended feature?')
+                 cal_catch = np.flipud(cal_catch)
+                 
 
-       
-       cal_catch[cal_catch<1] = -9999
+       cal_catch[cal_catch < 1] = -9999
        cal_catch[np.isnan(cal_catch) == 1] = -9999
        wr_id = np.unique(wr)
-       
+      
        id_error_wr = []
        cal_catch_error_wr = []
        output_message = []

Diff for: src/lisfloodutilities/writers/pcr.py

@@ -139,15 +139,17 @@ def _get_geo_transform_from_coords(self):
         # 4 - 0.0
         # 5 - n-s pixel resolution (negative value)
         #
-        w_e_resolution = np.round(self._coordinates['x'][1] - self._coordinates['x'][0], 3)
-        n_s_resolution = np.round(self._coordinates['y'][1] - self._coordinates['y'][0], 3)
+        n_s_resolution =(self._coordinates['y'][-1] - self._coordinates['y'][0])/(len(self._coordinates['y'][:])-1.0)
+        w_e_resolution =(self._coordinates['x'][-1] - self._coordinates['x'][0])/(len(self._coordinates['x'][:])-1.0)
         # considering coordinates at center of the cell, we must subtract half step for top left points
         if (self._coordinates['y'][0]<self._coordinates['y'][-1]):
            top=self._coordinates['y'][-1]
            ns=-1
         else:
            top=self._coordinates['y'][0]
            ns=1        
-        top_left_x = np.round(self._coordinates['x'][0], 2) - w_e_resolution / 2
-        top_left_y = np.round(top, 2) - ns*n_s_resolution / 2
-        return top_left_x, w_e_resolution, 0.0, top_left_y, 0.0, n_s_resolution
+        top_left_x = self._coordinates['x'][0] - w_e_resolution / 2 
+        top_left_y = top - ns*n_s_resolution / 2 
+        return top_left_x, w_e_resolution, 0.0, top_left_y, 0.0, n_s_resolution        
+        
+