Merge pull request #32 from swiss-territorial-data-lab/ch/add_FP_EPT

Ch/add fp ept

Author: clmarmy

examples/mineral-extract-sites-detection/config_det.yaml

@@ -19,7 +19,7 @@ generate_tilesets.py:
   datasets:
     aoi_tiles: output_det/tiles.geojson
     image_source:
-      type: XYZ     # supported values: 1. MIL = Map Image Layer 2. WMS 3. XYZ
+      type: XYZ     # supported values: 1. MIL = Map Image Layer 2. WMS 3. XYZ 4. FOLDER
       location: https://wmts.geo.admin.ch/1.0.0/ch.swisstopo.swissimage-product/default/2020/3857/{z}/{x}/{y}.jpeg
   output_folder: output_det
   tile_size: 256    # per side, in pixels

examples/mineral-extract-sites-detection/config_trne.yaml

@@ -2,25 +2,34 @@
 prepare_data.py:  
   srs: EPSG:2056
   datasets:
-    shapefile: ./data/labels/tlm-hr-trn-topo.shp
+    shapefile: ./data/labels/tlm-hr-trn-topo.shp                   # GT labels
+    FP_shapefile: ./data/FP/FP_list.gpkg                           # FP labels
+    empty_tiles_aoi: ./data/AoI/AoI_2020.shp                       # AOI in which additional empty tiles can be selected. Only one 'empty_tiles' option can be selected  
+    # empty_tiles_shp: .data/EPT/<SHPFILE>                         # Provided shpfile of selected empty tiles. Only one 'empty_tiles' option can be selected  
   output_folder: ./output/output_trne
-  zoom_level: 16 
+  zoom_level: 16
 
 # 2-Fetch of tiles and split into 3 datasets: train, test, validation
 generate_tilesets.py:
   debug_mode: 
     enable: False  # sample of tiles
-    nb_tiles_max: 500
+    nb_tiles_max: 1000
   working_directory: output
   datasets:
     aoi_tiles: output_trne/tiles.geojson
     ground_truth_labels: output_trne/labels.geojson
+    FP_labels: output_trne/FP.geojson               
     image_source:
-      type: XYZ                                   # supported values: 1. MIL = Map Image Layer 2. WMS 3. XYZ 4. FOLDER
-      year: 2020                            # supported values: 1. multi-year (tiles of different year), 2. <year> (i.e. 2020)
+      type: XYZ           # supported values: 1. MIL = Map Image Layer 2. WMS 3. XYZ 4. FOLDER
+      year: 2020
       location: https://wmts.geo.admin.ch/1.0.0/ch.swisstopo.swissimage-product/default/{year}/3857/{z}/{x}/{y}.jpeg
+  empty_tiles:            # add empty tiles to datasets
+    tiles_frac: 0.5       # fraction (relative to the number of tiles intersecting labels) of empty tiles to add
+    frac_trn: 0.75        # fraction of empty tiles to add to the trn dataset, then the remaining tiles will be split in 2 and added to tst and val datasets
+    keep_oth_tiles: False # keep tiles in oth dataset not intersecting oth labels
   output_folder: output_trne
-  tile_size: 256      # per side, in pixels
+  tile_size: 256          # per side, in pixels
+  seed: 42
   overwrite: True
   n_jobs: 10
   COCO_metadata:
@@ -57,7 +66,7 @@ make_detections.py:
     tst: COCO_tst.json
   detectron2_config_file: ../../detectron2_config_dqry.yaml # path relative to the working_folder
   model_weights:
-    pth_file: ./logs/model_final.pth # trained model minimising the validation loss curve, monitor the training process via tensorboard (tensorboard --logdir </logs>)
+    pth_file: ./logs/model_0000999.pth # trained model minimising the validation loss curve, monitor the training process via tensorboard (tensorboard --logdir </logs>)
   image_metadata_json: img_metadata.json
   rdp_simplification:   # rdp = Ramer-Douglas-Peucker
     enabled: true

examples/mineral-extract-sites-detection/data/FP/FP_list.gpkg

@@ -11,6 +11,7 @@
 import geopandas as gpd
 import morecantile
 import pandas as pd
+from shapely.geometry import Polygon
 
 sys.path.insert(0, '../..')
 from helpers.misc import format_logger
@@ -21,6 +22,14 @@
 
 
 def add_tile_id(row):
+    """Attribute tile id
+
+    Args:
+        row (DataFrame): row of a given df
+
+    Returns:
+        DataFrame: row with addition 'id' column
+    """
 
     re_search = re.search('(x=(?P<x>\d*), y=(?P<y>\d*), z=(?P<z>\d*))', row.title)
     if 'year' in row.keys():
@@ -31,15 +40,60 @@ def add_tile_id(row):
     return row
 
 
+def aoi_tiling(gdf):
+    """Tiling of an AoI
+
+    Args:
+        gdf (GeoDataFrame): gdf containing all the bbox boundary coordinates
+
+    Returns:
+        Geodataframe: gdf containing the tiles shape of the bbox of the AoI
+    """
+
+    # Grid definition
+    tms = morecantile.tms.get('WebMercatorQuad')    # epsg:3857
+
+    tiles_all = [] 
+    for boundary in gdf.itertuples():
+        coords = (boundary.minx, boundary.miny, boundary.maxx, boundary.maxy)      
+        tiles = gpd.GeoDataFrame.from_features([tms.feature(x, projected=False) for x in tms.tiles(*coords, zooms=[ZOOM_LEVEL])]) 
+        tiles.set_crs(epsg=4326, inplace=True)
+        tiles_all.append(tiles)
+    tiles_all_gdf = gpd.GeoDataFrame(pd.concat(tiles_all, ignore_index=True))
+
+    return tiles_all_gdf
+
+
+def bbox(bounds):
+    """Get a vector bounding box of a 2D shape
+
+    Args:
+        bounds (array): minx, miny, maxx, maxy of the bounding box
+
+    Returns:
+        geometry (Polygon): polygon geometry of the bounding box
+    """
+
+    minx = bounds[0]
+    miny = bounds[1]
+    maxx = bounds[2]
+    maxy = bounds[3]
+
+    return Polygon([[minx, miny],
+                    [maxx, miny],
+                    [maxx, maxy],
+                    [minx, maxy]])
+
+
 if __name__ == "__main__":
 
     # Start chronometer
     tic = time.time()
     logger.info('Starting...')
 
     # Argument and parameter specification
-    parser = argparse.ArgumentParser(description='The script prepares the Mineral Extraction Sites dataset to be processed by the object-detector scripts')
-    parser.add_argument('config_file', type=str, help='Framework configuration file')
+    parser = argparse.ArgumentParser(description="The script prepares the Mineral Extraction Sites dataset to be processed by the object-detector scripts")
+    parser.add_argument('config_file', type=str, help="Framework configuration file")
     args = parser.parse_args()
 
     logger.info(f"Using {args.config_file} as config file.")
@@ -50,6 +104,19 @@ def add_tile_id(row):
     # Load input parameters
     OUTPUT_DIR = cfg['output_folder']
     SHPFILE = cfg['datasets']['shapefile']
+    FP_SHPFILE = cfg['datasets']['FP_shapefile'] if 'FP_shapefile' in cfg['datasets'].keys() else None
+    if 'empty_tiles_aoi' in cfg['datasets'].keys() and 'empty_tiles_shp' in cfg['datasets'].keys():
+        logger.error("Choose one option between providing an AoI shapefile ('empty_tiles_aoi') in which empty tiles will be selected and a shapefile with selected empty tiles ('empty_tiles_shp')")
+        sys.exit(1)    
+    if 'empty_tiles_aoi' in cfg['datasets'].keys():
+        EPT_SHPFILE = cfg['datasets']['empty_tiles_aoi']
+        EPT = 'aoi'
+    elif 'empty_tiles_shp' in cfg['datasets'].keys():
+        EPT_SHPFILE = cfg['datasets']['empty_tiles_shp'] 
+        EPT = 'shp'
+    else:
+        EPT_SHPFILE = None
+        EPT = None
     ZOOM_LEVEL = cfg['zoom_level']
 
     # Create an output directory in case it doesn't exist
@@ -60,70 +127,124 @@ def add_tile_id(row):
 
     # Prepare the tiles
 
-    ## Convert datasets shapefiles into geojson format
-    logger.info('Convert labels shapefile into GeoJSON format (EPSG:4326)...')
+    # Convert datasets shapefiles into geojson format
+    logger.info("Convert labels shapefile into GeoJSON format (EPSG:4326)...")
     labels = gpd.read_file(SHPFILE)
     labels_4326 = labels.to_crs(epsg=4326)
     labels_4326['CATEGORY'] = 'quarry'
     labels_4326['SUPERCATEGORY'] = 'land usage'
 
     nb_labels = len(labels)
-    logger.info(f'There is/are {nb_labels} polygons in {SHPFILE}')
+    logger.info(f"There are {nb_labels} polygons in {SHPFILE}")
 
-    label_filename = 'labels.geojson'
-    label_filepath = os.path.join(OUTPUT_DIR, label_filename)
-    labels_4326.to_file(label_filepath, driver='GeoJSON')
-    written_files.append(label_filepath)  
-    logger.success(f"{DONE_MSG} A file was written: {label_filepath}")
+    filename = 'labels.geojson'
+    filepath = os.path.join(OUTPUT_DIR, filename)
+    labels_4326.to_file(filepath, driver='GeoJSON')
+    written_files.append(filepath)  
+    logger.success(f"{DONE_MSG} A file was written: {filepath}")
 
-    logger.info('Creating tiles for the Area of Interest (AoI)...')   
-    
-    # Grid definition
-    tms = morecantile.tms.get('WebMercatorQuad')    # epsg:3857
+    # Add FP labels if it exists
+    if FP_SHPFILE:
+        fp_labels = gpd.read_file(FP_SHPFILE)
+        fp_labels_4326 = fp_labels.to_crs(epsg=4326)
+        fp_labels_4326['CATEGORY'] = 'quarry'
+        fp_labels_4326['SUPERCATEGORY'] = 'land usage'
+
+        nb_fp_labels = len(fp_labels)
+        logger.info(f"There are {nb_fp_labels} polygons in {FP_SHPFILE}")
+
+        filename = 'FP.geojson'
+        filepath = os.path.join(OUTPUT_DIR, filename)
+        fp_labels_4326.to_file(filepath, driver='GeoJSON')
+        written_files.append(filepath)  
+        logger.success(f"{DONE_MSG} A file was written: {filepath}")
+
+        labels_4326 = pd.concat([labels_4326, fp_labels_4326], ignore_index=True)
+    else:
+        labels_4326 = labels_4326
+
+    # Keep only label boundary geometry info (minx, miny, maxx, maxy) 
+    logger.info("- Get the label boundaries")  
+    boundaries_df = labels_4326.bounds
+
+    # Get the boundaries for all the labels (minx, miny, maxx, maxy) 
+    global_boundaries_gdf = labels_4326.dissolve() if len(labels_4326) > 0 else labels_4326
+    labels_bbox = bbox(global_boundaries_gdf.iloc[0].geometry.bounds)
 
-    # New gpd with only labels geometric info (minx, miny, maxx, maxy) 
-    logger.info('- Get geometric boundaries of the labels')  
-    label_boundaries_df = labels_4326.bounds
-
-    # Iterate on geometric coordinates to defined tiles for a given label at a given zoom level
-    # A gpd is created for each label and are then concatenate into a single gpd 
-    logger.info('- Compute tiles for each labels geometry') 
-    tiles_4326_all = [] 
-
-    for label_boundary in label_boundaries_df.itertuples():
-        coords = (label_boundary.minx, label_boundary.miny, label_boundary.maxx, label_boundary.maxy)   
-        tiles_4326 = gpd.GeoDataFrame.from_features([tms.feature(x, projected=False) for x in tms.tiles(*coords, zooms=[ZOOM_LEVEL])])   
-        tiles_4326.set_crs(epsg=4326, inplace=True)
-        tiles_4326_all.append(tiles_4326)
-    tiles_4326_aoi = gpd.GeoDataFrame(pd.concat(tiles_4326_all, ignore_index=True))
-
-    # Remove unrelevant tiles and reorganised the data set:
-    logger.info('- Remove duplicated tiles and tiles that are not intersecting labels') 
-
-    # - Keep only tiles that are actually intersecting labels
-    labels_4326.rename(columns={'FID': 'id_aoi'}, inplace=True)
-    tiles_4326 = gpd.sjoin(tiles_4326_aoi, labels_4326, how='inner')
-
-    # - Remove duplicated tiles
-    if nb_labels > 1:
-        tiles_4326.drop_duplicates(['title', 'year'] if 'year' in tiles_4326.keys() else 'title', inplace=True)
-
-    # - Remove useless columns, reset feature id and redefine it according to xyz format  
-    logger.info('- Add tile IDs and reorganise data set')
-    tiles_4326 = tiles_4326[['geometry', 'title', 'year'] if 'year' in tiles_4326.keys() else ['geometry', 'title']].copy()
-    tiles_4326.reset_index(drop=True, inplace=True)
+    # Get tiles for a given AoI from which empty tiles will be selected when the images are retrieved
+    # Get tiles for a given AoI from which empty tiles will be selected when the images are retrieved
+    if EPT_SHPFILE:
+        EPT_aoi = gpd.read_file(EPT_SHPFILE)
+        EPT_aoi_4326 = EPT_aoi.to_crs(epsg=4326)
+        
+        if EPT == 'aoi':
+            logger.info("- Get AoI boundaries")  
+            EPT_aoi_boundaries_df = EPT_aoi_4326.bounds
+
+            # Get the boundaries for all the AoI (minx, miny, maxx, maxy) 
+            EPT_aoi_boundaries_gdf = EPT_aoi_4326.dissolve() if len(EPT_aoi_4326) > 0 else EPT_aoi_4326
+            aoi_bbox = bbox(EPT_aoi_boundaries_gdf.iloc[0].geometry.bounds)
+            aoi_bbox_contains = aoi_bbox.contains(labels_bbox)
+
+            if aoi_bbox_contains:
+                logger.info("- The surface area occupied by the bbox of the AoI used to find empty tiles is bigger than the label's one. The AoI boundaries will be used for tiling") 
+                boundaries_df = EPT_aoi_boundaries_df
+            else:
+                logger.info("- The surface area occupied by the bbox of the AoI used to find empty tiles is smaller than the label's one. Both the AoI and labels area will be used for tiling") 
+                # Get tiles coordinates and shapes
+                empty_tiles_4326_all = aoi_tiling(EPT_aoi_boundaries_df)
+                # Delete tiles outside of the AoI limits 
+                empty_tiles_4326_aoi = gpd.sjoin(empty_tiles_4326_all, EPT_aoi_4326, how='inner', lsuffix='ept_tiles', rsuffix='ept_aoi')
+        elif EPT == 'shp':
+            empty_tiles_4326_aoi = EPT_aoi_4326
+            aoi_bbox = None
+            aoi_bbox_contains = False
+
+    logger.info("Creating tiles for the Area of Interest (AoI)...")   
+
+    # Get tiles coordinates and shapes
+    tiles_4326_aoi = aoi_tiling(boundaries_df)
+
+    if EPT_SHPFILE and aoi_bbox_contains:
+        # Delete tiles outside of the AoI limits 
+        tiles_4326_aoi = gpd.sjoin(tiles_4326_aoi, EPT_aoi_4326, how='inner', lsuffix='ept_tiles', rsuffix='ept_aoi', predicate='intersects')
+
+    # Compute labels intersecting tiles 
+    tiles_gt_4326 = gpd.sjoin(tiles_4326_aoi, labels_4326, how='inner', predicate='intersects')
+    tiles_gt_4326.drop_duplicates('title', inplace=True)
+    logger.info(f"- Number of tiles intersecting GT labels = {len(tiles_gt_4326)}")
+    if FP_SHPFILE:
+        tiles_fp_4326 = gpd.sjoin(tiles_4326_aoi, fp_labels_4326, how='inner', predicate='intersects')
+        tiles_fp_4326.drop_duplicates('title', inplace=True)
+        logger.info(f"- Number of tiles intersecting FP labels = {len(tiles_fp_4326)}")
+
+    if not EPT_SHPFILE or EPT_SHPFILE and aoi_bbox_contains == False:
+        # Keep only tiles intersecting labels 
+        tiles_4326_aoi = gpd.sjoin(tiles_4326_aoi, labels_4326, how='inner', predicate='intersects')
+        tiles_4326_aoi.drop_duplicates('title', inplace=True)
+ 
+    # Get all the tiles in one gdf 
+    if EPT_SHPFILE and aoi_bbox_contains == False:
+        logger.info("- Add label tiles to empty AoI tiles") 
+        tiles_4326_all = pd.concat([tiles_4326_aoi, empty_tiles_4326_aoi])
+    else: 
+        tiles_4326_all = tiles_4326_aoi
 
-    # Add the ID column
-    tiles_4326 = tiles_4326.apply(add_tile_id, axis=1)
+    tiles_4326_all.drop_duplicates('title', inplace=True)
+
+    # Add tile IDs and reorganise data set
+    tiles_4326_all = tiles_4326_all[['geometry', 'title']].copy()
+    tiles_4326_all.reset_index(drop=True, inplace=True)
+    tiles_4326_all = tiles_4326_all.apply(add_tile_id, axis=1)
 
-    nb_tiles = len(tiles_4326)
-    logger.info(f'There was/were {nb_tiles} tiles created')
+    nb_tiles = len(tiles_4326_all)
+    logger.info(f"There were {nb_tiles} tiles created")
 
-    # Export tiles to GeoJSON
-    logger.info('Export tiles to GeoJSON (EPSG:4326)...')  
+    # Save tile shapefile
+    logger.info("Export tiles to GeoJSON (EPSG:4326)...")  
     tile_filename = 'tiles.geojson'
     tile_filepath = os.path.join(OUTPUT_DIR, tile_filename)
-    tiles_4326.to_file(tile_filepath, driver='GeoJSON')
+    tiles_4326_all.to_file(tile_filepath, driver='GeoJSON')
     written_files.append(tile_filepath)  
     logger.success(f"{DONE_MSG} A file was written: {tile_filepath}")
 

examples/mineral-extract-sites-detection/prepare_data.py

@@ -31,6 +31,10 @@ generate_tilesets.py:
       year: 2018
       location: https://wmts.geo.admin.ch/1.0.0/ch.swisstopo.swissimage-product/default/{year}/3857/{z}/{x}/{y}.jpeg
       srs: "EPSG:3857"
+  # empty_tiles:            # add empty tiles to datasets
+  #   tiles_frac: 0.5       # fraction (relative to the number of tiles intersecting labels) of empty tiles to add
+  #   frac_trn: 0.75        # fraction of empty tiles to add to the trn dataset, then the remaining tiles will be split in 2 and added to tst and val datasets
+  #   keep_oth_tiles: False # keep tiles in oth dataset not intersecting oth labels
   output_folder: .
   tile_size: 256 # per side, in pixels
   overwrite: False

examples/road-surface-classification/config_rs.yaml

@@ -9,16 +9,20 @@ prepare_data.py:
 generate_tilesets.py:
   debug_mode: 
     enable: False  # sample of tiles
-    nb_tiles_max: 100
+    nb_tiles_max: 1000
   working_directory: .
   datasets:
     aoi_tiles: output_GE/aoi_z18_tiles.geojson
     ground_truth_labels: output_GE/ground_truth_labels.geojson
     other_labels: output_GE/other_labels.geojson
     image_source:
-      type: MIL # supported values: 1. MIL = Map Image Layer 2. WMS
+      type: MIL # supported values: 1. MIL = Map Image Layer 2. WMS 3. XYZ 4. FOLDER
       location: https://raster.sitg.ge.ch/arcgis/rest/services/ORTHOPHOTOS_2018_EPSG2056/MapServer
       srs: "EPSG:3857"
+  # empty_tiles:            # add empty tiles to datasets
+  #   tiles_frac: 0.5       # fraction (relative to the number of tiles intersecting labels) of empty tiles to add
+  #   frac_trn: 0.75        # fraction of empty tiles to add to the trn dataset, then the remaining tiles will be split in 2 and added to tst and val datasets
+  #   keep_oth_tiles: True  # keep tiles in oth dataset not intersecting oth labels
   output_folder: output_GE
   tile_size: 256 # per side, in pixels
   overwrite: False

examples/swimming-pool-detection/GE/config_GE.yaml

@@ -12,7 +12,7 @@
 
 from loguru import logger
 
-sys.path.insert(0, '.')
+sys.path.insert(1, '../../..')
 from helpers.misc import format_logger
 
 logger = format_logger(logger)