Update "lfcoords" to include the reservoirs flag

bin/lfcoords

@@ -8,7 +8,7 @@ src_dir = os.path.normpath(os.path.join(current_dir, '../src/'))
 if os.path.exists(src_dir):
     sys.path.append(src_dir)
 
-from lisfloodutilities.lfcoords import main_script
+from lisfloodutilities.lfcoords.lfcoords import main_script
 
 if __name__ == '__main__':
     main_script()

src/lisfloodutilities/lfcoords/__init__.py

@@ -19,7 +19,7 @@ def __init__(self, config_file):
             config = yaml.load(ymlfile, Loader=yaml.FullLoader)
 
         # input
-        self.STATIONS = Path(config['input']['stations'])
+        self.POINTS = Path(config['input']['points'])
         self.LDD_FINE = Path(config['input']['ldd_fine'])
         self.UPSTREAM_FINE = Path(config['input']['upstream_fine'])
         self.LDD_COARSE = Path(config['input']['ldd_coarse'])

src/lisfloodutilities/lfcoords/coarser_grid.py

@@ -13,7 +13,7 @@
 warnings.filterwarnings("ignore")
 
 from lisfloodutilities.lfcoords import Config
-from lisfloodutilities.lfcoords.utils import catchment_polygon
+from lisfloodutilities.lfcoords.utils import catchment_polygon, downstream_pixel
 
 # set logger
 logging.basicConfig(level=logging.INFO,
@@ -22,28 +22,31 @@
 
 def coordinates_coarse(
     cfg: Config,
-    stations: pd.DataFrame,
+    points: pd.DataFrame,
+    reservoirs: bool = False,
     save: bool = True
 ) -> Optional[pd.DataFrame]:
     """
-    Transforms station coordinates from a high-resolution grid to a corresponding location in a coarser grid, aiming to match the shape of the catchment area derived from the high-resolution map. It updates the station coordinates and exports the catchment areas as shapefiles in the coarser grid.
+    Transforms point coordinates from a high-resolution grid to a corresponding location in a coarser grid, aiming to match the shape of the catchment area derived from the high-resolution map. It updates the point coordinates and exports the catchment areas as shapefiles in the coarser grid.
 
     The function reads the upstream area map and local drainage direction (LDD) map in the coarse grid. It then finds the pixel in the coarse grid that best matches the catchment shape derived from the high-resolution map. The match is evaluated based on the intersection-over-union of catchment shapes and the ratio of upstream areas.
 
     Parameters:
     -----------
     cfg: Config
         Configuration object containing file paths and parameters specified in the configuration file.
-    stations: pandas.DataFrame
-        DataFrame containing station coordinates and upstream areas in the finer grid. It's the result of coarse_grid.coarse_grid()
+    points: pandas.DataFrame
+        DataFrame containing point coordinates and upstream areas in the finer grid. It's the result of coarse_grid.coarse_grid()
+    reservoirs: bool
+        Whether the points are reservoirs or not. If True, the resulting coordinates refer to one pixel downstream of the actual solution, a deviation required by the LISFLOOD reservoir simulation
     save: bool
-        If True, the updated stations table is saved to a CSV file.
-        If False, the updated stations DataFrame is returned without saving.
+        If True, the updated points table is saved to a CSV file.
+        If False, the updated points DataFrame is returned without saving.
 
     Returns:
     --------
-    staions: pandas.DataFrame
-        If save is False, returns a pandas DataFrame with updated station coordinates and upstream areas in the coarser grid. Otherwise, the function returns None, and the results are saved directly to a CSV file.
+    points: pandas.DataFrame
+        If save is False, returns a pandas DataFrame with updated point coordinates and upstream areas in the coarser grid. Otherwise, the function returns None, and the results are saved directly to a CSV file.
     """
 
     ### READ INPUTS
@@ -56,6 +59,9 @@ def coordinates_coarse(
     ldd_coarse = rioxarray.open_rasterio(cfg.LDD_COARSE).squeeze(dim='band')
     logger.info(f'Map of local drainage directions correctly read: {cfg.LDD_COARSE}')
 
+    # copy of the points dataframe
+    points_ = points.copy()
+
 
     ### PROCESSING
 
@@ -75,13 +81,13 @@ def coordinates_coarse(
     suffix_coarse = f'{cellsize_arcmin}min'
     logger.info(f'Coarse resolution is {cellsize_arcmin} arcminutes')
 
-    # extract resolution of the finer grid from 'stations'
-    suffix_fine = sorted(set([col.split('_')[1] for col in stations.columns if '_' in col]))[0]
+    # extract resolution of the finer grid from 'points'
+    suffix_fine = sorted(set([col.split('_')[1] for col in points.columns if '_' in col]))[0]
     cols_fine = [f'{col}_{suffix_fine}' for col in ['area', 'lat', 'lon']]
 
-    # add new columns to 'stations'
+    # add new columns to 'points'
     cols_coarse = [f'{col}_{suffix_coarse}' for col in ['area', 'lat', 'lon']]
-    stations[cols_coarse] = np.nan
+    points_[cols_coarse] = np.nan
 
     # output folders
     SHAPE_FOLDER_FINE = cfg.SHAPE_FOLDER / suffix_fine
@@ -90,7 +96,7 @@ def coordinates_coarse(
 
     # search range of 5x5 array -> this is where the best point can be found in the coarse grid
     rangexy = np.linspace(-2, 2, 5) * cellsize # arcmin
-    for ID, attrs in tqdm(stations.iterrows(), total=stations.shape[0], desc='stations'):
+    for ID, attrs in tqdm(points_.iterrows(), total=points_.shape[0], desc='points'):
 
         # real upstream area
         area_ref = attrs['area']
@@ -99,7 +105,7 @@ def coordinates_coarse(
         lat_fine, lon_fine, area_fine = attrs[[f'{col}_{suffix_fine}' for col in ['lat', 'lon', 'area']]]
 
         if (area_ref < cfg.MIN_AREA) or (area_fine < cfg.MIN_AREA):
-            logger.warning(f'The catchment area of station {ID} is smaller than the minimum of {cfg.MIN_AREA} km2')
+            logger.warning(f'The catchment area of point {ID} is smaller than the minimum of {cfg.MIN_AREA} km2')
             continue
 
         # import shapefile of catchment polygon
@@ -183,14 +189,18 @@ def coordinates_coarse(
         basin_coarse.to_file(output_shp)
         logger.info(f'Catchment {ID} exported as shapefile: {output_shp}')
 
-        # update new columns in 'stations'
-        stations.loc[ID, cols_coarse] = [int(area_coarse), round(lat_coarse, 6), round(lon_coarse, 6)]
+        # move the result one pixel downstream, in case of reservoir
+        if reservoirs:
+            lat_coarse, lon_coarse = downstream_pixel(lat_coarse, lon_coarse, upstream_coarse)
+
+        # update new columns in 'points_'
+        points_.loc[ID, cols_coarse] = [int(area_coarse), round(lat_coarse, 6), round(lon_coarse, 6)]
 
     # return/save
-    stations.sort_index(axis=1, inplace=True)
+    points_.sort_index(axis=1, inplace=True)
     if save:
-        output_csv = cfg.STATIONS.parent / f'{cfg.STATIONS.stem}_{suffix_coarse}.csv'
-        stations.to_csv(output_csv)
-        logger.info(f'The updated stations table in the coarser grid has been exported to: {output_csv}')
+        output_csv = cfg.POINTS.parent / f'{cfg.POINTS.stem}_{suffix_coarse}.csv'
+        points_.to_csv(output_csv)
+        logger.info(f'The updated points table in the coarser grid has been exported to: {output_csv}')
     else:
-        return stations
+        return points_

src/lisfloodutilities/lfcoords/finer_grid.py

@@ -22,28 +22,28 @@ def coordinates_fine(
     save: bool = True
 ) -> Optional[pd.DataFrame]:
     """
-    Processes station coordinates to find the most accurate pixel in a high-resolution map, 
-    based on a reference value of catchment area. It updates the station coordinates and 
+    Processes point coordinates to find the most accurate pixel in a high-resolution map, 
+    based on a reference value of catchment area. It updates the point coordinates and 
     exports the catchment areas as shapefiles.
 
     The function reads the upstream area map and local drainage direction (LDD) map in 
-    fine resolution, as well as the station coordinates with their reference upstream areas.
+    fine resolution, as well as the point coordinates with their reference upstream areas.
     It then iteratively searches for the best-matching pixel within an increasing search 
     range and applies penalties and factors to determine the closest match. The function 
-    creates a boolean map of each station's catchment area and vectorizes it into a 
+    creates a boolean map of each point's catchment area and vectorizes it into a 
     GeoDataFrame for export as a shapefile.
 
     Parameters:
     -----------
     cfg: Config
         Configuration object containing file paths and parameters specified in  the configuration file.
     save: bool
-        If True, the updated stations table is saved to a CSV file. If False, the updated stations DataFrame is returned without saving.
+        If True, the updated points table is saved to a CSV file. If False, the updated points DataFrame is returned without saving.
         
     Returns:
     --------
-    stations: pandas.DataFrame
-        If save is False, returns a pandas DataFrame with updated station coordinates and upstream areas in the finer grid. Otherwise, the function returns None, and the results are saved directly to a CSV file.
+    points: pandas.DataFrame
+        If save is False, returns a pandas DataFrame with updated point coordinates and upstream areas in the finer grid. Otherwise, the function returns None, and the results are saved directly to a CSV file.
     """
 
     # READ INPUTS
@@ -56,9 +56,9 @@ def coordinates_fine(
     ldd_fine = rioxarray.open_rasterio(cfg.LDD_FINE).squeeze(dim='band')
     logger.info(f'Map of local drainage directions corretly read: {cfg.LDD_FINE}')
 
-    # read stations text file
-    stations = pd.read_csv(cfg.STATIONS, index_col='ID')
-    logger.info(f'Table of stations correctly read: {cfg.STATIONS}')
+    # read points text file
+    points = pd.read_csv(cfg.POINTS, index_col='ID')
+    logger.info(f'Table of points correctly read: {cfg.POINTS}')
 
 
     # PROCESSING
@@ -69,9 +69,9 @@ def coordinates_fine(
     suffix_fine = f'{cellsize_arcsec}sec'
     logger.info(f'The resolution of the finer grid is {cellsize_arcsec} arcseconds')
 
-    # add columns to the table of stations
-    new_cols = sorted([f'{col}_{suffix_fine}' for col in stations.columns])
-    stations[new_cols] = np.nan
+    # add columns to the table of points
+    new_cols = sorted([f'{col}_{suffix_fine}' for col in points.columns])
+    points[new_cols] = np.nan
 
     # create river network
     fdir_fine = pyflwdir.from_array(ldd_fine.data,
@@ -84,7 +84,7 @@ def coordinates_fine(
     SHAPE_FOLDER_FINE = cfg.SHAPE_FOLDER / suffix_fine
     SHAPE_FOLDER_FINE.mkdir(parents=True, exist_ok=True)
 
-    for ID, attrs in tqdm(stations.iterrows(), total=stations.shape[0], desc='stations'):  
+    for ID, attrs in tqdm(points.iterrows(), total=points.shape[0], desc='points'):  
 
         # reference coordinates and upstream area
         lat_ref, lon_ref, area_ref = attrs[['lat', 'lon', 'area']]
@@ -100,8 +100,8 @@ def coordinates_fine(
             if error <= max_error:
                 break
 
-        # update new columns in 'stations'
-        stations.loc[ID, new_cols] = [int(upstream_fine.sel(y=lat, x=lon).item()), round(lat, 6), round(lon, 6)]
+        # update new columns in 'points'
+        points.loc[ID, new_cols] = [int(upstream_fine.sel(y=lat, x=lon).item()), round(lat, 6), round(lon, 6)]
 
         # boolean map of the catchment associated to the corrected coordinates
         basin_arr = fdir_fine.basins(xy=(lon, lat)).astype(np.int32)
@@ -121,10 +121,10 @@ def coordinates_fine(
         logger.info(f'Catchment {ID} exported as shapefile: {output_file}')
 
     # return/save
-    stations.sort_index(axis=1, inplace=True)
+    points.sort_index(axis=1, inplace=True)
     if save:
-        output_csv = cfg.STATIONS.parent / f'{cfg.STATIONS.stem}_{suffix_fine}.csv'
-        stations.to_csv(output_csv)
-        logger.info(f'The updated stations table in the finer grid has been exported to: {output_csv}')
+        output_csv = cfg.POINTS.parent / f'{cfg.POINTS.stem}_{suffix_fine}.csv'
+        points.to_csv(output_csv)
+        logger.info(f'The updated points table in the finer grid has been exported to: {output_csv}')
     else: 
-        return stations
+        return points

src/lisfloodutilities/lfcoords/lfcoords.py

@@ -17,14 +17,14 @@
 import logging
 
 from lisfloodutilities.lfcoords import Config
-from lilisfloodutilities.lfcoords.finer_grid import coordinates_fine
-from llisfloodutilities.lfcoords.coarser_grid import coordinates_coarse
+from lisfloodutilities.lfcoords.finer_grid import coordinates_fine
+from lisfloodutilities.lfcoords.coarser_grid import coordinates_coarse
 
 def main(argv=sys.argv):
     prog = os.path.basename(argv[0])
     parser = argparse.ArgumentParser(
         description="""
-        Correct the coordinates of a set of stations to match the river network in the
+        Correct the coordinates of a set of points to match the river network in the
         LISFLOOD static map.
         First, it uses a reference value of catchment area to find the most accurate
         pixel in a high-resolution map.
@@ -33,11 +33,14 @@ def main(argv=sys.argv):
         """,
         prog=prog
     )
-    parser.add_argument('-c', '--config-file', type=str, required=True, help='Path to the YML configuration file')
+    parser.add_argument('-c', '--config-file', type=str, required=True,
+                        help='Path to the YML configuration file')
+    parser.add_argument('-r', '--reservoirs', action='store_true', default=False,
+                        help='The input points are reservoirs')
     args = parser.parse_args()
 
     # create logger
-    logger = logging.getLogger('correct-coordinates')
+    logger = logging.getLogger('lfcoords')
     logger.setLevel(logging.INFO)
     logger.propagate = False
     log_format = logging.Formatter('%(asctime)s | %(levelname)s | %(name)s | %(message)s', datefmt='%Y-%m-%d %H:%M:%S')
@@ -55,20 +58,20 @@ def main(argv=sys.argv):
 
     # find coordinates in high resolution
     try:
-        stations_HR = coordinates_fine(cfg, save=False)
+        points_HR = coordinates_fine(cfg, save=False)
     except Exception as e:
         logger.error(f'Locating the points in the finer grid: {e}')
         sys.exit(2)
 
     # find coordinates in LISFLOOD
     try:
-        coordinates_coarse(cfg, stations_HR, save=True)
+        coordinates_coarse(cfg, points_HR, reservoirs=args.reservoirs, save=True)
     except Exception as e:
         logger.error(f'Locating the points in the finer grid: {e}')
         sys.exit(3)
 
 def main_script():
-    sys.exist(main())
+    sys.exit(main())
 
 if __name__ == "__main__":
     main_script()

src/lisfloodutilities/lfcoords/utils.py

@@ -87,6 +87,50 @@ def find_pixel(
 
 
 
+def downstream_pixel(
+    lat: float,
+    lon: float,
+    upArea: xr.DataArray
+) -> (float, float):
+    """It finds the downstream coordinates of a given point
+    
+    Parameteres:
+    ------------
+    lat: float
+        latitude of the input point
+    lon: float
+        longitued of the input point
+    upArea: xarray.DataArray
+        map of upstream area
+        
+    Returns:
+    --------
+    lat: float
+        latitude of the inmediate downstream pixel
+    lon: float
+        longitued of the inmediate downstream pixel
+    """
+
+    # upstream area of the input coordinates
+    area = upArea.sel(x=lon, y=lat, method='nearest').item()
+
+    # spatial resolution of the input map
+    resolution = np.mean(np.diff(upArea.x.values))
+
+    # window around the input pixel
+    window = np.array([-1.5 * resolution, 1.5 * resolution])
+    upArea_ = upArea.sel(y=slice(*window[::-1] + lat)).sel(x=slice(*window + lon))
+
+    # remove pixels with area equal or smaller than the input pixel
+    mask = upArea_.where(upArea_ > area, np.nan)
+
+    # from the remaining, find pixel with the smallest upstream area
+    pixel = upArea_.where(upArea_ == mask.min(), drop=True)
+
+    return round(pixel.y.item(), 6), round(pixel.x.item(), 6)
+
+
+
 def catchment_polygon(
     data: np.ndarray,
     transform: Affine,

tests/data/lfcoords/config.yml

@@ -1,5 +1,5 @@
 input:
-    stations: data/lfcoords/points.csv
+    points: data/lfcoords/points.csv
     ldd_fine: data/lfcoords/MERIT/ldd_3sec.tif
     upstream_fine: data/lfcoords/MERIT/uparea_3sec.tif
     ldd_coarse: data/lfcoords/EFAS/ldd_1min.nc