Merge pull request #10 from IGNF/debug/indices-map

Debug/indices map

Author: leavauchier

CHANGELOG.md

@@ -1,5 +1,7 @@
 # CHANGELOG
 
+- génération de la carte d'indice même quand il n'y a pas de points à ajouter
+
 ## 1.1.1
 - lint
 - ajout de pre-commit hooks pour appliquer le lint au moment des commits

environment.yml

@@ -9,6 +9,9 @@ dependencies:
   - numpy
   - geopandas==0.*
   - shapely>=2.0.3
+  - gdal
+  - pdal>=2.8
+  - python-pdal
   - rasterio
   # ------------- logging ------------- #
   - loguru
@@ -18,4 +21,6 @@ dependencies:
   # ----------- linting --------------- #
   - pre-commit
   - black
-  - flake8
+  - flake8
+  - pip:
+    - ign-pdal-tools==1.8.1

patchwork/indices_map.py

@@ -8,38 +8,36 @@
 from rasterio.transform import from_origin
 
 from patchwork.constants import PATCH_X_STR, PATCH_Y_STR
-from patchwork.tools import get_tile_origin_from_pointcloud
 
 
-def create_indices_grid(config: DictConfig, df_points: DataFrame) -> np.ndarray:
+def create_indices_grid(config: DictConfig, df_points: DataFrame, corner_x: int, corner_y: int) -> np.ndarray:
     """create a binary grid matching the tile the points of df_points are from, where each patch is equal to:
     1 if the patch has at least one point of df_points
     0 if the patch has no point from df_points
     """
     size_grid = int(config.TILE_SIZE / config.PATCH_SIZE)
 
-    corner_x, corner_y = get_tile_origin_from_pointcloud(config, df_points)
+    grid = np.zeros((size_grid, size_grid))
 
-    list_coordinates_x = np.int32((df_points.x - corner_x) / config.PATCH_SIZE)
-    list_coordinates_y = np.int32((corner_y - df_points.y) / config.PATCH_SIZE)
+    if not df_points.empty:
+        list_coordinates_x = np.int32((df_points.x - corner_x) / config.PATCH_SIZE)
+        list_coordinates_y = np.int32((corner_y - df_points.y) / config.PATCH_SIZE)
 
-    # edge cases where points are exactly on the... edge of the tile, but still valid
-    list_coordinates_x[list_coordinates_x == size_grid] = size_grid - 1
-    list_coordinates_y[list_coordinates_y == size_grid] = size_grid - 1
+        # edge cases where points are exactly on the... edge of the tile, but still valid
+        list_coordinates_x[list_coordinates_x == size_grid] = size_grid - 1
+        list_coordinates_y[list_coordinates_y == size_grid] = size_grid - 1
 
-    grid = np.zeros((size_grid, size_grid))
+        grid[list_coordinates_x, list_coordinates_y] = 1
 
-    grid[list_coordinates_x, list_coordinates_y] = 1
     return grid.transpose()
 
 
-def create_indices_map(config: DictConfig, df_points: DataFrame):
+def create_indices_map(config: DictConfig, df_points: DataFrame, corner_x: int, corner_y: int):
     """
     Save a binary grid for the tile into a geotiff
     """
-    corner_x, corner_y = get_tile_origin_from_pointcloud(config, df_points)
 
-    grid = create_indices_grid(config, df_points)
+    grid = create_indices_grid(config, df_points, corner_x, corner_y)
     os.makedirs(config.filepath.OUTPUT_INDICES_MAP_DIR, exist_ok=True)
     output_indices_map_path = os.path.join(
         config.filepath.OUTPUT_INDICES_MAP_DIR, config.filepath.OUTPUT_INDICES_MAP_NAME

patchwork/patchwork.py

@@ -8,6 +8,7 @@
 import pandas as pd
 from laspy import LasReader, ScaleAwarePointRecord
 from omegaconf import DictConfig
+from pdaltools.las_info import get_tile_origin_using_header_info
 
 import patchwork.constants as c
 from patchwork.indices_map import create_indices_map
@@ -125,20 +126,20 @@ def get_field_from_header(las_file: LasReader) -> List[str]:
     return [dimension.name.lower() for dimension in header.point_format.dimensions]
 
 
-def test_field_exists(file_path: str, colmun: str) -> bool:
+def test_field_exists(file_path: str, column: str) -> bool:
     output_file = laspy.read(file_path)
-    return colmun in get_field_from_header(output_file)
+    return column in get_field_from_header(output_file)
 
 
 def append_points(config: DictConfig, extra_points: pd.DataFrame):
     # get field to copy :
     recipient_filepath = os.path.join(config.filepath.RECIPIENT_DIRECTORY, config.filepath.RECIPIENT_NAME)
-    ouput_filepath = os.path.join(config.filepath.OUTPUT_DIR, config.filepath.OUTPUT_NAME)
+    output_filepath = os.path.join(config.filepath.OUTPUT_DIR, config.filepath.OUTPUT_NAME)
     with laspy.open(recipient_filepath) as recipient_file:
         recipient_fields_list = get_field_from_header(recipient_file)
 
     # get fields that are in the donor file we can transmit to the recipient without problem
-    # classification is in the fields to exclude because it will be copy in a special way
+    # classification is in the fields to exclude because it will be copied in a special way
     fields_to_exclude = [
         c.PATCH_X_STR,
         c.PATCH_Y_STR,
@@ -152,7 +153,7 @@ def append_points(config: DictConfig, extra_points: pd.DataFrame):
         if (field.lower() in extra_points.columns) and (field.lower() not in fields_to_exclude)
     ]
 
-    copy2(recipient_filepath, ouput_filepath)
+    copy2(recipient_filepath, output_filepath)
 
     if len(extra_points) == 0:  # if no point to add, the job is done after copying the recipient file
         return
@@ -165,11 +166,11 @@ def append_points(config: DictConfig, extra_points: pd.DataFrame):
                              column name in {recipient_filepath}"
             )
         new_column_type = get_type(config.NEW_COLUMN_SIZE)
-        output_las = laspy.read(ouput_filepath)
+        output_las = laspy.read(output_filepath)
         output_las.add_extra_dim(laspy.ExtraBytesParams(name=config.NEW_COLUMN, type=new_column_type))
-        output_las.write(ouput_filepath)
+        output_las.write(output_filepath)
 
-    with laspy.open(ouput_filepath, mode="a") as output_las:
+    with laspy.open(output_filepath, mode="a") as output_las:
         # put in a new table all extra points and their values on the fields we want to keep
         new_points = laspy.ScaleAwarePointRecord.zeros(extra_points.shape[0], header=output_las.header)
         for field in fields_to_keep:
@@ -227,12 +228,15 @@ def get_donor_path(config: DictConfig) -> Tuple[str, str]:
 
 def patchwork(config: DictConfig):
     _, donor_name = get_donor_path(config)
-    if not donor_name:  # if no matching donor, we simply copy the recipient to the output without doing anything
-        recipient_filepath = os.path.join(config.filepath.RECIPIENT_DIRECTORY, config.filepath.RECIPIENT_NAME)
-        ouput_filepath = os.path.join(config.filepath.OUTPUT_DIR, config.filepath.OUTPUT_NAME)
-        copy2(recipient_filepath, ouput_filepath)
-        return
+    recipient_filepath = os.path.join(config.filepath.RECIPIENT_DIRECTORY, config.filepath.RECIPIENT_NAME)
+    if donor_name:
+        complementary_bd_points = get_complementary_points(config)
+        append_points(config, complementary_bd_points)
+
+    else:  # if no matching donor, we simply copy the recipient to the output without doing anything
+        output_filepath = os.path.join(config.filepath.OUTPUT_DIR, config.filepath.OUTPUT_NAME)
+        copy2(recipient_filepath, output_filepath)
+        complementary_bd_points = pd.DataFrame()  # No points to add
 
-    complementary_bd_points = get_complementary_points(config)
-    append_points(config, complementary_bd_points)
-    create_indices_map(config, complementary_bd_points)
+    corner_x, corner_y = get_tile_origin_using_header_info(filename=recipient_filepath, tile_width=config.TILE_SIZE)
+    create_indices_map(config, complementary_bd_points, corner_x, corner_y)

test/test_indices_map.py

@@ -15,6 +15,8 @@
 
 PATCH_SIZE = 1
 TILE_SIZE = 3
+CORNER_X = 0
+CORNER_Y = 3
 
 DATA_POINTS = {"x": [0.0, 1.5, 3, 1.5, 2.5], "y": [0.0, 0.5, 0.5, 1.5, 3]}
 # we want y=0 at the bottom, but in a ndarray it's at the top, so grid['y'] = SIZE_Y - data_points['y']
@@ -28,7 +30,7 @@ def test_create_indices_points():
             config_name="configs_patchwork.yaml", overrides=[f"PATCH_SIZE={PATCH_SIZE}", f"TILE_SIZE={TILE_SIZE}"]
         )
         df_points = pd.DataFrame(data=DATA_POINTS)
-        grid = create_indices_grid(config, df_points)
+        grid = create_indices_grid(config, df_points, CORNER_X, CORNER_Y)
 
         grid = grid.transpose()  # indices aren't read the way we want otherwise
 
@@ -54,18 +56,42 @@ def test_create_indices_map(tmp_path_factory):
         )
 
         df_points = pd.DataFrame(data=DATA_POINTS)
-        create_indices_map(config, df_points)
+        create_indices_map(config, df_points, CORNER_X, CORNER_Y)
         raster = rs.open(os.path.join(tmp_file_dir, tmp_file_name))
         grid = raster.read()
 
         grid = grid.transpose()  # indices aren't read the way we want otherwise
+        print(grid)
 
         for point in POINTS_IN_GRID:
             assert grid[point] == 1
         for point in POINTS_NOT_IN_GRID:
             assert grid[point] == 0
 
 
+def test_create_indices_map_no_added_points(tmp_path_factory):
+    tmp_file_dir = tmp_path_factory.mktemp("data")
+    tmp_file_name = "empty_indices.tif"
+
+    with initialize(version_base="1.2", config_path="../configs"):
+        config = compose(
+            config_name="configs_patchwork.yaml",
+            overrides=[
+                f"PATCH_SIZE={PATCH_SIZE}",
+                f"TILE_SIZE={TILE_SIZE}",
+                f"filepath.OUTPUT_INDICES_MAP_DIR={tmp_file_dir}",
+                f"filepath.OUTPUT_INDICES_MAP_NAME={tmp_file_name}",
+            ],
+        )
+
+        df_points = pd.DataFrame()  # Empty dataframe (no points to add)
+        create_indices_map(config, df_points, CORNER_X, CORNER_Y)
+        raster = rs.open(os.path.join(tmp_file_dir, tmp_file_name))
+        grid = raster.read()
+
+        assert np.all(grid == 0)
+
+
 def test_read_indices_map(tmp_path_factory):
     tmp_file_dir = tmp_path_factory.mktemp("data")
     tmp_file_name = "indices.tif"

test/test_patchwork.py

@@ -16,6 +16,7 @@
     get_field_from_header,
     get_selected_classes_points,
     get_type,
+    patchwork,
 )
 from patchwork.tools import get_tile_origin_from_pointcloud
 
@@ -366,3 +367,66 @@ def test_get_donor_path(tmp_path_factory):
         donor_dir, donor_name = get_donor_path(config)
         assert donor_dir == DONOR_MORE_FIELDS_TEST_DIR
         assert donor_name == DONOR_MORE_FIELDS_TEST_NAME
+
+
+def test_patchwork_default(tmp_path_factory):
+    tmp_file_dir = tmp_path_factory.mktemp("data")
+    tmp_output_las_name = "result_patchwork.laz"
+    tmp_output_indices_map_name = "result_patchwerk_indices.tif"
+
+    with initialize(version_base="1.2", config_path="../configs"):
+        config = compose(
+            config_name="configs_patchwork.yaml",
+            overrides=[
+                f"filepath.RECIPIENT_DIRECTORY={RECIPIENT_TEST_DIR}",
+                f"filepath.RECIPIENT_NAME={RECIPIENT_TEST_NAME}",
+                f"filepath.DONOR_DIRECTORY={DONOR_TEST_DIR}",
+                f"filepath.DONOR_NAME={DONOR_TEST_NAME}",
+                f"filepath.OUTPUT_DIR={tmp_file_dir}",
+                f"filepath.OUTPUT_NAME={tmp_output_las_name}",
+                f"filepath.OUTPUT_INDICES_MAP_DIR={tmp_file_dir}",
+                f"filepath.OUTPUT_INDICES_MAP_NAME={tmp_output_indices_map_name}",
+            ],
+        )
+        patchwork(config)
+        recipient_path = os.path.join(config.filepath.RECIPIENT_DIRECTORY, config.filepath.RECIPIENT_NAME)
+        output_path = os.path.join(tmp_file_dir, tmp_output_las_name)
+        indices_map_path = os.path.join(tmp_file_dir, tmp_output_indices_map_name)
+        assert os.path.isfile(output_path)
+        assert os.path.isfile(indices_map_path)
+        with laspy.open(recipient_path) as las_file:
+            recipient_points = las_file.read().points
+        with laspy.open(output_path) as las_file:
+            output_points = las_file.read().points
+        assert len(output_points) > len(recipient_points)
+
+
+def test_patchwork_empty_donor(tmp_path_factory):
+    tmp_file_dir = tmp_path_factory.mktemp("data")
+    tmp_output_las_name = "result_patchwork.laz"
+    tmp_output_indices_map_name = "result_patchwerk_indices.tif"
+
+    with initialize(version_base="1.2", config_path="../configs"):
+        config = compose(
+            config_name="configs_patchwork.yaml",
+            overrides=[
+                f"filepath.RECIPIENT_DIRECTORY={RECIPIENT_TEST_DIR}",
+                f"filepath.RECIPIENT_NAME={RECIPIENT_TEST_NAME}",
+                f"filepath.OUTPUT_DIR={tmp_file_dir}",
+                f"filepath.OUTPUT_NAME={tmp_output_las_name}",
+                f"filepath.OUTPUT_INDICES_MAP_DIR={tmp_file_dir}",
+                f"filepath.OUTPUT_INDICES_MAP_NAME={tmp_output_indices_map_name}",
+            ],
+        )
+        patchwork(config)
+        recipient_path = os.path.join(config.filepath.RECIPIENT_DIRECTORY, config.filepath.RECIPIENT_NAME)
+        output_path = os.path.join(tmp_file_dir, tmp_output_las_name)
+        indices_map_path = os.path.join(tmp_file_dir, tmp_output_indices_map_name)
+        assert os.path.isfile(output_path)
+        assert os.path.isfile(indices_map_path)
+        with laspy.open(recipient_path) as las_file:
+            recipient_points = las_file.read().points
+        with laspy.open(output_path) as las_file:
+            output_points = las_file.read().points
+
+        assert len(output_points) == len(recipient_points)