diff --git a/README.md b/README.md index 824ea29..5dfc652 100644 --- a/README.md +++ b/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') ``` diff --git a/src/lisfloodutilities/mctrivers/mctrivers.py b/src/lisfloodutilities/mctrivers/mctrivers.py index fa5079c..d6f4b2b 100644 --- a/src/lisfloodutilities/mctrivers/mctrivers.py +++ b/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 diff --git a/tests/data/mctrivers/LF_ETRS89_UseCase/reference/mctmask.nc b/tests/data/mctrivers/LF_ETRS89_UseCase/reference/mctmask.nc index 6cb0b4f..2f6c6e9 100644 Binary files a/tests/data/mctrivers/LF_ETRS89_UseCase/reference/mctmask.nc and b/tests/data/mctrivers/LF_ETRS89_UseCase/reference/mctmask.nc differ diff --git a/tests/data/mctrivers/LF_lat_lon_UseCase/reference/mctmask.nc b/tests/data/mctrivers/LF_lat_lon_UseCase/reference/mctmask.nc index c24f074..25b6b80 100644 Binary files a/tests/data/mctrivers/LF_lat_lon_UseCase/reference/mctmask.nc and b/tests/data/mctrivers/LF_lat_lon_UseCase/reference/mctmask.nc differ diff --git a/tests/test_mctrivers.py b/tests/test_mctrivers.py index 74d4449..0f9d9a8 100644 --- a/tests/test_mctrivers.py +++ b/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() \ No newline at end of file + self.run(0.001, 5, 500*10**6, ['lat', 'lon'], 'LF_lat_lon_UseCase') \ No newline at end of file