Merge pull request #69 from NikosMastrantonas/make_mct_rivers_mask

debug tests for mctrivers

README.md

@@ -776,16 +776,34 @@ The tool can take the following additional input arguments:
 - `-S`, `--slope`: Riverbed slope threshold to use MCT diffusive wave routing (default:  0.001)
 - `-N`, `--nloops`: Number of consecutive downstream grid cells that also need to comply with the slope requirement for including a grid cell in the MCT rivers mask (default: 5)
 - `-U`, `--minuparea`: Minimum upstream drainage area for a pixel to be included in the MCT rivers mask (uses the same units as in the -u file) (default: 0)
-- `-E`, `--coordsnames`: Coordinates names for lat, lon (in this order with space!) used in the netcdf files
+- `-E`, `--coordsnames`: Coordinates names for lat, lon (in this order with space!) used in the netcdf files. The function checks for 3 commonly used names (x, lon, rlon for longitudes, and y, lat, rlat for latitudes). Therefere, it is recommended to keep the default value.
 
 The tool generates the following outputs:
 
 - `-O`, `--outputfilename`: Output file containing the rivers mask where LISFLOOD can use the MCT diffusive wave routing  (default: chanmct.nc)
 
-Example of command that will generate an MCT rivers mask with pixels where riverbed slope < 0.001, drainage area > 500 kms and at least 5 downstream pixels meet the same two conditions, considering the units of the upArea.nc file are given in kms:
+It can be used either from command line, or also as a python function. Below follow some examples:
+
+Example of command that will generate an MCT rivers mask with pixels where riverbed slope < 0.001 (same as default), drainage area > 500 kms and at least 5 (same as default) downstream pixels meet the same two conditions, considering the units of the upArea.nc file are given in kms:
 
 ```bash
-mctrivers -i changrad.nc -l ldd.nc -u upArea.nc -O chanmct.nc [-m mask.nc -E y x -S 0.001 -N 5 -U 500 optional] 
+mctrivers -i changrad.nc -l ldd.nc -u upArea.nc -O chanmct.nc -U 500
+```
+
+```python
+from lisfloodutilities.mctrivers.mctrivers import mct_mask
+mct_mask(channels_slope_file='changrad.nc', ldd_file='ldd.nc', uparea_file='upArea.nc', minuparea=500, outputfile='chanmct.nc')
+```
+
+Example of command that will generate an MCT rivers mask with pixels where riverbed slope < 0.0005, drainage area > 0 (same as default) kms and at least 3 downstream pixels meet the same two conditions. Also a mask (mask.nc) will be used, and the coords names in the nc files are "Lat1", "Lon1" for lat, lon respectively:
+
+```bash
+mctrivers -i changrad.nc -l ldd.nc -u upArea.nc -O chanmct.nc -m mask.nc -E Lat1 Lon1 -S 0.0005 -N 3 
+```
+
+```python
+from lisfloodutilities.mctrivers.mctrivers import mct_mask
+mct_mask(channels_slope_file='changrad.nc', ldd_file='ldd.nc', uparea_file='upArea.nc', mask_file='mask.nc', slp_threshold=0.0005, nloops=3, coords=["Lat1", "Lon1"], outputfile='chanmct.nc')
 ```
 
 

src/lisfloodutilities/mctrivers/mctrivers.py

@@ -62,14 +62,14 @@ def mct_mask(channels_slope_file, ldd_file, uparea_file, mask_file='',
     slp_threshold: Riverbed slope threshold to use MCT diffusive wave routing (default: 0.001)
     nloops: Number of consecutive downstream grid cells that also need to comply with the slope requirement for including a grid cell in the MCT rivers mask (default: 5)
     minuparea: Minimum upstream drainage area for a pixel to be included in the MCT rivers mask (uses the same units as in the -u file) (default: 0)
-    coords_names: Coordinates names for lat, lon (in this order as list) used in the the netcdf files (default: 'None': checks for commonly used names)
+    coords_names: Coordinates names for lat, lon (in this order as list) used in the the netcdf files (default: 'None'; checks for commonly used names ['x', 'lon', 'rlon'], similar for lat names)
     outputfile: Output file containing the rivers mask where LISFLOOD can use the MCT diffusive wave routing (default: chanmct.nc)
     
     Example for generating an MCT rivers mask with pixels where riverbed slope < 0.001, drainage area > 500 kms and at least 5 downstream pixels meet the same 
     two conditions, considering the units of the upArea.nc file are given in kms:
     
     mct_mask(channels_slope_file='changrad.nc', ldd_file='ldd.nc', uparea_file='upArea.nc', mask_file='mask.nc', 
-             slp_threshold=0.001, nloops=5, minuparea=0, coords_names=['y' , 'x'], 
+             slp_threshold=0.001, nloops=5, minuparea=500, coords_names=['y' , 'x'], 
              outputfile='chanmct.nc')
     """
     # ---------------- Read LDD (Note that for EFAS5 there is small shift of values for CH)
@@ -78,7 +78,7 @@ def mct_mask(channels_slope_file, ldd_file, uparea_file, mask_file='',
     # ---------------- Auxiliary variables
     x_checks = ['lon', 'x', 'rlon']
     y_checks = ['lat', 'y', 'rlat']
-    if coords_names[0] == "None":
+    if coords_names == "None":
         x_proj = set(list(LD.coords)) & set(x_checks)
         y_proj = set(list(LD.coords)) & set(y_checks)
 
@@ -167,7 +167,7 @@ def mct_mask(channels_slope_file, ldd_file, uparea_file, mask_file='',
         MX = MX.fillna(0)*0+1
 
     # use the exact same coords from channel slope file, just in case there are precision differences
-    MX = MX.assign_coords(x_proj=x_all, y_proj=y_all)
+    MX = MX.assign_coords({x_proj: x_all, y_proj: y_all})
     LD.close()  # close ther LD file, after the check of mask availability
 
     # ---------------- Loop on the basin pixels to find how many MCT pixels they have downstream

tests/test_mctrivers.py

@@ -2,61 +2,64 @@
 import shutil
 import xarray as xr
 
-from lisfloodutilities.compare.nc import NetCDFComparator
-
 from lisfloodutilities.mctrivers.mctrivers import mct_mask
 
 
 def mk_path_out(p):
+    'make folder for storing output data'
     path_out = os.path.join(os.path.dirname(__file__), p)
     if os.path.exists(path_out):
-        shutil.rmtree(path_out)
+        shutil.rmtree(path_out, ignore_errors=True)
     os.mkdir(path_out)
-    return path_out
+
 
 class TestMctMask():
 
     case_dir = os.path.join(os.path.dirname(__file__), 'data', 'mctrivers')
 
-    def run(self, slp_threshold, nloops, minuparea, coords_names, type):
-
+    def run(self, slp_threshold, nloops, minuparea, coords_names, case_name):
+        'generate the MCT mask'
         # setting
-        self.out_path_ref = os.path.join(self.case_dir, type, 'reference')
-        self.out_path_run = os.path.join(self.case_dir, type, 'out')
+        self.out_path_ref = os.path.join(self.case_dir, case_name, 'reference')
+        self.out_path_run = os.path.join(self.case_dir, case_name, 'out')
         mk_path_out(self.out_path_run)
 
-        channels_slope_file = os.path.join(self.case_dir, type, 'changrad.nc')
-        ldd_file = os.path.join(self.case_dir, type, 'ldd.nc')
-        uparea_file = os.path.join(self.case_dir, type, 'upArea.nc')
-        mask_file = os.path.join(self.case_dir, type, 'mask.nc')
-        outputfile = os.path.join(self.out_path_run, 'chanmct.nc')
+        channels_slope_file = os.path.join(self.case_dir, case_name, 'changrad.nc')
+        ldd_file = os.path.join(self.case_dir, case_name, 'ldd.nc')
+        uparea_file = os.path.join(self.case_dir, case_name, 'upArea.nc')
+        mask_file = os.path.join(self.case_dir, case_name, 'mask.nc')
+        outputfile = os.path.join(self.out_path_run, 'mctmask.nc')
 
         # generate the mct river mask
-        mct_mask(channels_slope_file, ldd_file, uparea_file, mask_file, slp_threshold, nloops, minuparea,outputfile)
-
-        # compare results with reference
-        nc_comparator = NetCDFComparator(mask_file, array_equal=True)
-        nc_comparator.compare_dirs(self.out_path_run, self.out_path_ref)
-
-    def teardown_method(self):
-        print('Cleaning directories')
-        out_path = os.path.join(self.case_dir, type, 'out')
-        if os.path.exists(out_path) and os.path.isdir(out_path):
-            shutil.rmtree(out_path, ignore_errors=True)
+        mct_mask(channels_slope_file, ldd_file, uparea_file, mask_file, slp_threshold, nloops, minuparea, coords_names, outputfile)
+
+        # compare the generated mask with the reference one
+        ref_file = self.out_path_ref+'/mctmask.nc'
+        reference = xr.open_dataset(ref_file)
+        out_file = self.out_path_run+'/mctmask.nc'
+        generated = xr.open_dataset(out_file)
+        # check if same based on https://docs.xarray.dev/en/stable/generated/xarray.DataArray.equals.html
+        all_equal = reference.equals(generated)
+        generated.close()  # needs to be closed otherwise the out folder can't be deleted
+
+        if not all_equal:
+            fail_message = f'Test for mct river mask generation for {case_name} failed. Please check differences between' 
+            fail_message += f' the generated mask "{out_file}" and the expected mask "{ref_file}".'
+            assert all_equal, fail_message
+        else:
+            # if equal just delete the out folder
+            shutil.rmtree(self.out_path_run, ignore_errors=True)
 
 
 class TestMctrivers(TestMctMask):
 
     # slp_threshold = 0.001
     # nloops = 5
-    # minuparea = 0
+    # minuparea = 500*10**6
     # coords_names = 'None'
 
     def test_mctrivers_etrs89(self):
-        self.run(0.001, 5, 0, 'None', 'LF_ETRS89_UseCase')
+        self.run(0.001, 5, 500*10**6, 'None', 'LF_ETRS89_UseCase')
 
     def test_mctrivers_latlon(self):
-        self.run( 0.001, 5, 0, 'x' 'y', 'LF_lat_lon_UseCase')
-
-    def cleaning(self,):
-        self.teardown_method()
+        self.run(0.001, 5, 500*10**6, ['lat', 'lon'], 'LF_lat_lon_UseCase')