Merge pull request #70 from NikosMastrantonas/make_mct_rivers_mask

Make mct rivers mask

Author: doc78

.gitignore

@@ -1,4 +1,5 @@
 # Created by .ignore support plugin (hsz.mobi)
+.ipynb*
 .idea
 *.nc
 **/__pycache__/
@@ -57,4 +58,4 @@ tests/data/gridding/meteo_out/**/*.tiff
 .project
 .pydevproject
 
-.ipynb_checkpoints/
+.ipynb_checkpoints/

README.md

@@ -778,7 +778,7 @@ The tool can take the following additional input arguments:
 - `-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. 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:
+The tool generates the following outputs (when called from command line as main script):
 
 - `-O`, `--outputfilename`: Output file containing the rivers mask where LISFLOOD can use the MCT diffusive wave routing  (default: chanmct.nc)
 
@@ -791,8 +791,9 @@ mctrivers -i changrad.nc -l ldd.nc -u upArea.nc -O chanmct.nc -U 500
 ```
 
 ```python
+# no data are saved when called inside 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')
+mct_mask_ds = mct_mask(channels_slope_file='changrad.nc', ldd_file='ldd.nc', uparea_file='upArea.nc', minuparea=500)
 ```
 
 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:
@@ -803,7 +804,7 @@ mctrivers -i changrad.nc -l ldd.nc -u upArea.nc -O chanmct.nc -m mask.nc -E Lat1
 
 ```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')
+mct_mask_ds = 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"])
 ```
 
 

src/lisfloodutilities/mctrivers/mctrivers.py

@@ -43,8 +43,7 @@ def getarg():
 
 
 def mct_mask(channels_slope_file, ldd_file, uparea_file, mask_file='', 
-             slp_threshold=0.001, nloops=5, minuparea=0, coords_names='None', 
-             outputfile='chanmct.nc'):
+             slp_threshold=0.001, nloops=5, minuparea=0, coords_names='None'):
     """
     
     Builds a mask of mild sloping rivers for use in LISFLOOD with MCT diffusive river routing. It takes LISFLOOD channels slope map (changrad.nc), the LDD (ldd.nc), 
@@ -57,7 +56,7 @@ def mct_mask(channels_slope_file, ldd_file, uparea_file, mask_file='',
     
     channels_slope_file: LISFLOOD channels gradient map (changrad.nc)
     ldd_file: LISFLOOD local drain direction file (ldd.nc)
-    uparea_file: LISFLOOD Uustream area file (upArea.nc)
+    uparea_file: LISFLOOD Upstream area file (upArea.nc)
     mask_file: a mask nc file; if not given (default) all cells are considered valid.
     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)
@@ -69,8 +68,7 @@ def mct_mask(channels_slope_file, ldd_file, uparea_file, mask_file='',
     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=500, coords_names=['y' , 'x'], 
-             outputfile='chanmct.nc')
+             slp_threshold=0.001, nloops=5, minuparea=500, coords_names=['y' , 'x'])
     """
     # ---------------- Read LDD (Note that for EFAS5 there is small shift of values for CH)
     LD = xr.open_dataset(ldd_file)
@@ -98,7 +96,7 @@ def mct_mask(channels_slope_file, ldd_file, uparea_file, mask_file='',
     # proprocess CH dataset for having correct format
     CH = xr.open_dataset(channels_slope_file)
     old_name = [i for i in list(CH.data_vars) if sorted(CH[i].dims)==sorted([x_proj, y_proj])]
-    CH = CH.rename({old_name[0]: "changrad"})  # only 1 variable complies with above check
+    CH = CH.rename({old_name[0]: "changrad"})  # only 1 variable complies with the above check
     CH['changrad'] = CH['changrad'].transpose(y_proj, x_proj)  # make sure dims order is as pcraster needs
 
     # ---------------- Set clone map for pcraster
@@ -149,7 +147,7 @@ def mct_mask(channels_slope_file, ldd_file, uparea_file, mask_file='',
     # ---------------- Read upstream area
     UA = xr.open_dataset(uparea_file)
     old_name = [i for i in list(UA.data_vars) if sorted(UA[i].dims)==sorted([x_proj, y_proj])]
-    UA = UA.rename({old_name[0]: "domain"})['domain']  # only 1 variable complies with above if
+    UA = UA.rename({old_name[0]: "domain"})['domain']  # only 1 variable complies with the above if statement
 
     # convert the xarray to pcraster
     UA = UA>=minuparea # check that the area is over the minimum
@@ -160,11 +158,10 @@ def mct_mask(channels_slope_file, ldd_file, uparea_file, mask_file='',
     try:
         MX = xr.open_dataset(mask_file)
         old_name = [i for i in list(MX.data_vars) if sorted(MX[i].dims)==sorted([x_proj, y_proj])]
-        MX = MX.rename({old_name[0]: "domain"})['domain']  # only 1 variable complies with above if
+        MX = MX.rename({old_name[0]: "domain"})['domain']  # only 1 variable complies with the above if statement
     except:
         print(f'The given mask path {mask_file} is not a valid path. All domain read from LDD file {ldd_file} is considered vaid.')
         MX = LD.copy(deep=True)
-        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})
@@ -180,11 +177,19 @@ def mct_mask(channels_slope_file, ldd_file, uparea_file, mask_file='',
     # Loop nloops times and use downstream function to find out if each cell has nloop MCT cells downstream
     # Downstream function gives the value in the downstream pixel in a map:
     # here it gives 1 if the downstream pixel is Muskingum, zero otherwise.
+    # Note that for pits it always gives its own value so we need to mask out all pits as ann intermidiate step
+    
+    # modify data, so that the most downstream point is masked out (otherwise the below loop gives for all points the same value)
+    downstream_actual_mask_pcr = pcr.downstreamdist(ldd_pcr)
+    downstream_actual_mask_pcr = pcr.ifthenelse(downstream_actual_mask_pcr == 0, pcr.boolean(0), pcr.boolean(1))
+    downstream_actual_mask_pcr = pcr.scalar(downstream_actual_mask_pcr)
+    
     # The loop is used to count how many pixels are MCT downstream, as at each loop we move the values 1 pixel upstream
     # At the end of the loop, each element of the array has the number of downstream MCT pixels for that pixel
     for loops in range(0, nloops):
         # get the value on the downstream cell and put it in a mask
         downstream_cells_pcr = pcr.downstream(ldd_pcr, downstream_cells_pcr)
+        downstream_cells_pcr = downstream_cells_pcr*downstream_actual_mask_pcr
         sum_rivers_pcr = sum_rivers_pcr + downstream_cells_pcr
 
     # ---------------- Generate a new MCT rivers mask
@@ -207,16 +212,14 @@ def mct_mask(channels_slope_file, ldd_file, uparea_file, mask_file='',
     MCT.name = 'mct_mask'
 
     # mask final data with the mask_file
-    MCT = MCT.where(MX==1)
+    MCT = MCT.where(MX.notnull())
 
-    # lisflood does not read NaNs so the data are saved as boolean 0-1, with 0 being flagged as NaN for python reading
-    MCT.to_netcdf(outputfile, encoding={"mct_mask": {'_FillValue': 0, 'dtype': 'int8'}})
     return MCT
 
 
 def main():
-    'function for runnign from command line'
-    # ---------------- Read settings
+    'function for running from command line'
+    # ---------------- Read settingss
     args = getarg()
     channels_slope_file_arg = args.changradfile
     ldd_file_arg = args.LDDfile
@@ -228,10 +231,11 @@ def main():
     coords_names_arg = args.coordsnames
     outputfile_arg = args.outputfilename
 
-    mct_mask(channels_slope_file=channels_slope_file_arg, ldd_file=ldd_file_arg, uparea_file=uparea_file_arg, mask_file=mask_file_arg, 
-             slp_threshold=slp_threshold_arg, nloops=nloops_arg, minuparea=minuparea_arg, coords_names=coords_names_arg, 
-             outputfile=outputfile_arg)
-
+    mct_final = mct_mask(channels_slope_file=channels_slope_file_arg, ldd_file=ldd_file_arg, uparea_file=uparea_file_arg, mask_file=mask_file_arg, 
+                         slp_threshold=slp_threshold_arg, nloops=nloops_arg, minuparea=minuparea_arg, coords_names=coords_names_arg)
+    # lisflood does not read NaNs so the data are saved as boolean 0-1, with -1 being flagged as NaN for python reading
+    mct_final.to_netcdf(outputfile_arg, encoding={"mct_mask": {'_FillValue': -1, 'dtype': 'int8'}})
+    
 
 if __name__ == "__main__":
     main()

tests/data/mctrivers/LF_ETRS89_UseCase/reference/mctmask.nc

@@ -31,8 +31,9 @@ def run(self, slp_threshold, nloops, minuparea, coords_names, case_name):
         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, coords_names, outputfile)
-
+        mct_final = mct_mask(channels_slope_file, ldd_file, uparea_file, mask_file, slp_threshold, nloops, minuparea, coords_names)
+        mct_final.to_netcdf(outputfile, encoding={"mct_mask": {'_FillValue': 0, 'dtype': 'int8'}})
+        
         # compare the generated mask with the reference one
         ref_file = self.out_path_ref+'/mctmask.nc'
         reference = xr.open_dataset(ref_file)