Update "lfcoords" to create log file and check points before processing

src/lisfloodutilities/lfcoords/__init__.py

@@ -12,10 +12,7 @@
 import logging
 
 # set logger
-logging.basicConfig(level=logging.INFO,
-                    format='%(asctime)s | %(levelname)s | %(name)s | %(message)s', datefmt='%Y-%m-%d %H:%M:%S')
-logger = logging.getLogger(__name__)
-
+logger = logging.getLogger('lfcoords')
 
 
 class Config:
@@ -59,28 +56,27 @@ def update_config(
         fine_grid: xr.DataArray,
         coarse_grid: xr.DataArray
     ):
-    """It extracts the resolution of the finer and coarser grid, updates the respective attributes in the configuration object, and it creates the necessary structure of directories
-    
-    Parameters:
-    -----------
-    fine_grid: xarray.DataArray
-        Any map in the fine grid
-    coarse_grid: xarray.DataArray
-        Any map in the coarse grid    
-    """
-
-    # resolution of the finer grid
-    cellsize = np.mean(np.diff(fine_grid.x)) # degrees
-    cellsize_arcsec = int(np.round(cellsize * 3600, 0)) # arcsec
-    logger.info(f'The resolution of the finer grid is {cellsize_arcsec} arcseconds')
-    self.FINE_RESOLUTION = f'{cellsize_arcsec}sec'
-
-    # resolution of the input maps
-    cellsize = np.round(np.mean(np.diff(coarse_grid.x)), 6) # degrees
-    cellsize_arcmin = int(np.round(cellsize * 60, 0)) # arcmin
-    logger.info(f'The resolution of the coarser grid is {cellsize_arcmin} arcminutes')
-    self.COARSE_RESOLUTION = f'{cellsize_arcmin}min'
-
+        """It extracts the resolution of the finer and coarser grid, updates the respective attributes in the configuration object, and it creates the necessary structure of directories
+
+        Parameters:
+        -----------
+        fine_grid: xarray.DataArray
+            Any map in the fine grid
+        coarse_grid: xarray.DataArray
+            Any map in the coarse grid    
+        """
+
+        # resolution of the finer grid
+        cellsize = np.mean(np.diff(fine_grid.x)) # degrees
+        cellsize_arcsec = int(np.round(cellsize * 3600, 0)) # arcsec
+        logger.info(f'The resolution of the finer grid is {cellsize_arcsec} arcseconds')
+        self.FINE_RESOLUTION = f'{cellsize_arcsec}sec'
+
+        # resolution of the input maps
+        cellsize = np.round(np.mean(np.diff(coarse_grid.x)), 6) # degrees
+        cellsize_arcmin = int(np.round(cellsize * 60, 0)) # arcmin
+        logger.info(f'The resolution of the coarser grid is {cellsize_arcmin} arcminutes')
+        self.COARSE_RESOLUTION = f'{cellsize_arcmin}min'
 
 
 def read_input_files(
@@ -123,6 +119,8 @@ def read_input_files(
     points = pd.read_csv(cfg.POINTS, index_col='ID')
     points.columns = points.columns.str.lower()
     logger.info(f'Table of points correctly read: {cfg.POINTS}')
+    points = check_points(cfg, points, ldd_fine)
+
     # convert to geopandas and export as shapefile
     points = gpd.GeoDataFrame(points,
                               geometry=[Point(xy) for xy in zip(points['lon'], points['lat'])],
@@ -143,3 +141,49 @@ def read_input_files(
     cfg.update_config(ldd_fine, ldd_coarse)
 
     return inputs
+
+
+def check_points(
+    cfg: Config,
+    points: pd.DataFrame,
+    ldd: xr.DataArray    
+    ) -> pd.DataFrame:
+    """Removes input points with missing value, catchment area smaller than the predefined threshold, or outside the extent of the input map
+    
+    Parameters:
+    -----------
+    cfg: Config
+        Configuration object containing file paths and parameters specified in the configuration file.
+    points: pandas.DataFrame
+        Table of input points with fields 'lat', 'lon' and 'area' (km2)
+    ldd: xarray.DataArray
+        Map of local drainage directions
+        
+    Returns:
+    --------
+    points: pandas.DataFrame
+        The input table in which points with conflicts have been removed
+    """
+
+    # remove points with missing values
+    mask_nan = points.isnull().any(axis=1)
+    if mask_nan.sum() > 0:
+        points = points[~mask_nan]
+        logger.warning(f'{mask_nan.sum()} points were removed because of missing values')
+
+    # remove points with small catchment area
+    mask_area = points['area'] < cfg.MIN_AREA
+    if mask_area.sum() > 0:
+        points = points[~mask_area]
+        logger.info(f'{mask_area.sum()} points were removed due to their small catchment area')
+
+    # remove points outside the input LDD map
+    lon_min, lat_min, lon_max, lat_max = np.round(ldd.rio.bounds(), 6)
+    mask_lon = (points.lon < lon_min) | (points.lon > lon_max)
+    mask_lat = (points.lat < lat_min) | (points.lat > lat_max)
+    mask_extent = mask_lon | mask_lat
+    if mask_extent.sum() > 0:
+        points = points[~mask_extent]
+        logger.info(f'{mask_extent.sum()} points were removed because they are outside the input LDD map')
+
+    return points

src/lisfloodutilities/lfcoords/coarser_grid.py

@@ -18,9 +18,8 @@
 from lisfloodpreprocessing.utils import catchment_polygon, downstream_pixel
 
 # set logger
-logging.basicConfig(level=logging.INFO,
-                    format='%(asctime)s | %(levelname)s | %(name)s | %(message)s', datefmt='%Y-%m-%d %H:%M:%S')
-logger = logging.getLogger(__name__)
+logger = logging.getLogger('lfcoords')
+
 
 def coordinates_coarse(
     cfg: Config,
@@ -62,6 +61,7 @@ def coordinates_coarse(
     """
 
     points = points_fine.copy()
+    n_points = points.shape[0]
 
     # create river network
     fdir_coarse = pyflwdir.from_array(ldd_coarse.data,
@@ -71,7 +71,7 @@ def coordinates_coarse(
                                       latlon=True)
 
     # boundaries and resolution of the input maps
-    lon_min, lat_min, lon_max, lat_max = np.round(ldd_coarse.rio.bounds(), 6)
+    # lon_min, lat_min, lon_max, lat_max = np.round(ldd_coarse.rio.bounds(), 6)
     cellsize = np.round(np.mean(np.diff(ldd_coarse.x)), 6) # degrees
 
     # extract resolution of the finer grid from 'points'
@@ -84,16 +84,16 @@ 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
     polygons_coarse = []
-    for ID, attrs in tqdm(points.iterrows(), total=points.shape[0], desc='points'):
+    pbar = tqdm(points.iterrows(), total=n_points, desc='points')
+    for n, (ID, attrs) in enumerate(pbar, start=1):
 
         # real upstream area
         area_ref = attrs['area']
-
+            
         # coordinates and upstream area in the fine grid
         lat_fine, lon_fine, area_fine = attrs[[f'{col}_{cfg.FINE_RESOLUTION}' 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')
+        if area_fine < cfg.MIN_AREA:
+            logger.warning(f'Skipping point {ID} because its catchment area in the finer grid is smaller than {cfg.MIN_AREA} km2')
             continue
 
         # find ratio
@@ -170,6 +170,8 @@ def coordinates_coarse(
         # update new columns in 'points'
         points.loc[ID, cols_coarse] = [int(area_coarse), round(lat_coarse, 6), round(lon_coarse, 6)]
 
+        logger.info(f'Point {ID} ({n}/{n_points}) located in the coarser grid')
+
     # concatenate polygons shapefile
     polygons_coarse = pd.concat(polygons_coarse)
 

src/lisfloodutilities/lfcoords/finer_grid.py

@@ -18,9 +18,8 @@
 from lisfloodpreprocessing.utils import find_pixel, catchment_polygon 
 
 # set logger
-logging.basicConfig(level=logging.INFO,
-                    format='%(asctime)s | %(levelname)s | %(name)s | %(message)s', datefmt='%Y-%m-%d %H:%M:%S')
-logger = logging.getLogger(__name__)
+logger = logging.getLogger('lfcoords')
+
 
 def coordinates_fine(
     cfg: Config,
@@ -62,15 +61,12 @@ def coordinates_fine(
         A table with the catchment polygons in the finer grid.
     """
 
-    # resolution of the input map
-    cellsize = np.mean(np.diff(upstream_fine.x)) # degrees
-    cellsize_arcsec = int(np.round(cellsize * 3600, 0)) # arcsec
-    # cfg.FINE_RESOLUTION = f'{cellsize_arcsec}sec'
-    fine_resolution = f'{cellsize_arcsec}sec'
-
     # add columns to the table of points
-    new_cols = sorted([f'{col}_{fine_resolution}' for col in ['lat', 'lon', 'area']])
-    points[new_cols] = np.nan
+    points_fine = points.copy()
+    n_points = points.shape[0]
+    cols = ['lat', 'lon', 'area']
+    new_cols = sorted([f'{col}_{cfg.FINE_RESOLUTION}' for col in cols])
+    points_fine[new_cols] = np.nan
 
     # create river network
     fdir_fine = pyflwdir.from_array(ldd_fine.data,
@@ -79,11 +75,12 @@ def coordinates_fine(
                                     check_ftype=False,
                                     latlon=True)
 
-    polygons = []
-    for ID, attrs in tqdm(points.iterrows(), total=points.shape[0], desc='points'):  
+    polygons_fine = []
+    pbar = tqdm(points.iterrows(), total=points.shape[0], desc='points')
+    for n, (ID, attrs) in enumerate(pbar, start=1):  
 
         # reference coordinates and upstream area
-        lat_ref, lon_ref, area_ref = attrs[['lat', 'lon', 'area']]
+        lat_ref, lon_ref, area_ref = attrs[cols]
 
         # search new coordinates in an increasing range
         ranges = [55, 101, 151]
@@ -96,8 +93,8 @@ def coordinates_fine(
             if error <= max_error:
                 break
 
-        # 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)]
+        # update new columns in 'points_fine'
+        points_fine.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)
@@ -108,28 +105,30 @@ def coordinates_fine(
                                       crs=ldd_fine.rio.crs,
                                       name='ID')
         basin_gdf['ID'] = ID
+        basin_gdf[cols] = attrs[cols]
         basin_gdf.set_index('ID', inplace=True)
-        basin_gdf[attrs.index] = attrs.values
 
         # save polygon
-        polygons.append(basin_gdf)
+        polygons_fine.append(basin_gdf)
+
+        logger.info(f'Point {ID} ({n}/{n_points}) located in the finer grid')
 
     # concatenate polygons shapefile
-    polygons = pd.concat(polygons)
+    polygons_fine = pd.concat(polygons_fine)
 
     # convert points to geopandas
-    geometry = [Point(xy) for xy in zip(points[f'lon_{fine_resolution}'], points[f'lat_{fine_resolution}'])]
-    points = gpd.GeoDataFrame(points, geometry=geometry, crs=4326)
-    points.sort_index(axis=1, inplace=True)
+    geometry = [Point(xy) for xy in zip(points_fine[f'lon_{cfg.FINE_RESOLUTION}'], points_fine[f'lat_{cfg.FINE_RESOLUTION}'])]
+    points_fine = gpd.GeoDataFrame(points_fine, geometry=geometry, crs=4326)
+    points_fine.sort_index(axis=1, inplace=True)
 
     if save is True:
         # polygons
-        polygon_shp = cfg.OUTPUT_FOLDER / f'catchments_{fine_resolution}.shp'
-        polygons.to_file(polygon_shp)
+        polygon_shp = cfg.OUTPUT_FOLDER / f'catchments_{cfg.FINE_RESOLUTION}.shp'
+        polygons_fine.to_file(polygon_shp)
         logger.info(f'Catchments in the finer grid have been exported to: {polygon_shp}')
         # points
-        point_shp = cfg.OUTPUT_FOLDER / f'{cfg.POINTS.stem}_{fine_resolution}.shp'
-        points.to_file(point_shp)
+        point_shp = cfg.OUTPUT_FOLDER / f'{cfg.POINTS.stem}_{cfg.FINE_RESOLUTION}.shp'
+        points_fine.to_file(point_shp)
         logger.info(f'The updated points table in the finer grid has been exported to: {point_shp}')
 
-    return points, polygons
+    return points_fine, polygons_fine

src/lisfloodutilities/lfcoords/lfcoords.py

@@ -15,10 +15,12 @@
 import argparse
 import pandas as pd
 import logging
+from datetime import datetime
 
-from lisfloodutilities.lfcoords import Config
-from lisfloodutilities.lfcoords.finer_grid import coordinates_fine
-from lisfloodutilities.lfcoords.coarser_grid import coordinates_coarse
+from lisfloodpreprocessing import Config, read_input_files
+from lisfloodpreprocessing.utils import find_conflicts
+from lisfloodpreprocessing.finer_grid import coordinates_fine
+from lisfloodpreprocessing.coarser_grid import coordinates_coarse
 
 def main(argv=sys.argv):
     prog = os.path.basename(argv[0])
@@ -44,10 +46,16 @@ def main(argv=sys.argv):
     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')
+    # console handler
     c_handler = logging.StreamHandler()
     c_handler.setFormatter(log_format)
     c_handler.setLevel(logging.INFO)
     logger.addHandler(c_handler)
+    # File handler
+    f_handler = logging.FileHandler(f'lfcoords_{datetime.now():%Y%m%d%H%M}.log')
+    f_handler.setFormatter(log_format)
+    f_handler.setLevel(logging.INFO)
+    logger.addHandler(f_handler)
 
     # read configuration
     try: