Refactor occupancy map generation

coclico/metrics/commons.py

@@ -1,5 +1,8 @@
+import logging
 from typing import List, Tuple
 
+import numpy as np
+
 from coclico.config import composed_class_separator
 
 
@@ -48,3 +51,31 @@ def bounded_affine_function(coordinates_min: Tuple, coordinates_max: Tuple, x_qu
         y_query = a * x_query + b
 
     return y_query
+
+
+def get_raster_geometry_from_las_bounds(las_bounds: Tuple[float], pixel_size: float):
+    """Compute expected raster top left corner (cell center) and number of pixels from las min/max values and pixel
+    size
+
+    Args:
+        las_bounds (Tuple(float)): Las min/max values : (x_min, y_min, x_max, ymax)
+        pixel_size (_type_): Pixel size of the raster
+
+    Returns:
+        _type_: coordinates of the top-left corner, and number of pixels on each axis
+    """
+    x_min_las, y_min_las, x_max_las, y_max_las = las_bounds
+    x_min = np.round(x_min_las / pixel_size) * pixel_size
+    y_min = np.round(y_min_las / pixel_size) * pixel_size
+    x_max = np.round(x_max_las / pixel_size) * pixel_size
+    y_max = np.round(y_max_las / pixel_size) * pixel_size
+    logging.debug(
+        f"Raster min/max cell centers infered from in las file: min: ({x_min}, {y_min}) max:({x_max}, {y_max})"
+    )
+
+    # add 1 to x_pixel and y_pixel because the map should be able to include the extreme points
+    x_pixels = int(np.ceil((x_max - x_min) / pixel_size)) + 1
+    y_pixels = int(np.ceil((y_max - y_min) / pixel_size)) + 1
+    logging.debug(f"Raster number of pixels infered from in las file: ({x_pixels}, {y_pixels})")
+
+    return (x_min, y_max), (x_pixels, y_pixels)

coclico/metrics/occupancy_map.py

@@ -0,0 +1,107 @@
+import logging
+from typing import Tuple
+
+import laspy
+import numpy as np
+import rasterio
+
+from coclico.metrics.commons import (
+    get_raster_geometry_from_las_bounds,
+    split_composed_class,
+)
+
+
+def create_2d_binary_map(
+    xs: np.array,
+    ys: np.array,
+    pixel_size: float,
+    x_min: float,
+    y_max: float,
+    nb_pixels: Tuple[int, int] = (1000, 1000),
+) -> np.array:
+    """Create 2d binary occupancy map from points coordinates:
+    boolean 2d map with 1 where at least one point falls in the pixel, 0 everywhere else.
+
+    Args:
+        xs (np.array): vector of x coordinates of all points
+        ys (np.array): vector of y coordinates of all points
+        pixel_size (float): pixel size (in meters) of the output map
+        x_min (float): x coordinate (in meters) of the pixel center of the upper left corner of the map
+        y_max (float): y coordinate (in meters) of the pixel center of the upper left corner of the map
+        nb_pixels (float, optional): number of pixels on each axis in format (x, y). Defaults to (1000, 1000).
+
+    Returns:
+        np.array: Boolean output map
+    """
+    # numpy array is filled with (y, x) instead of (x, y)
+    grid = np.zeros((nb_pixels[1], nb_pixels[0]), dtype=bool)
+
+    for x, y in zip(xs, ys):
+        grid_x = min(int((x - (x_min - pixel_size / 2)) / pixel_size), nb_pixels[0] - 1)  # x_min is left pixel center
+        grid_y = min(int(((y_max + pixel_size / 2) - y) / pixel_size), nb_pixels[1] - 1)  # y_max is upper pixel center
+        grid[grid_y, grid_x] = True
+
+    return grid
+
+
+def create_multilayer_2d_occupancy_map(las_file, class_weights, output_tif, pixel_size):
+    """Create 2d occupancy map for each class that is in class_weights keys, and save result in a single output_tif
+    file with one layer per class (the classes are sorted alphabetically).
+
+    Args:
+        las_file (Path): path to the las file on which to generate occupancy map
+        class_weights (Dict): class weights dict (to know for which classes to generate the binary map)
+        output_tif (Path): path to output
+        pixel_size (float): size of the output raster pixels
+    """
+    with laspy.open(las_file) as f:
+        las = f.read()
+        xs, ys = las.x, las.y
+        classifs = las.classification
+        crs = las.header.parse_crs()
+
+    las_bounds = (np.min(xs), np.min(ys), np.max(xs), np.max(ys))
+
+    top_left, nb_pixels = get_raster_geometry_from_las_bounds(las_bounds, pixel_size)
+    x_min, y_max = top_left
+
+    def create_binary_map_from_class(class_key):
+        splitted_class_key = split_composed_class(class_key)
+        splitted_class_key_int = [int(ii) for ii in splitted_class_key]
+
+        map = create_2d_binary_map(
+            xs[np.isin(classifs, splitted_class_key_int)],
+            ys[np.isin(classifs, splitted_class_key_int)],
+            pixel_size,
+            x_min,
+            y_max,
+            nb_pixels,
+        )
+
+        return map
+
+    # get results for classes that are in weights dictionary (merged if necessary) in a 3d array
+    # which represents a raster with 1 layer per class
+    # keys are sorted to make sure that raster layers can be retrieved in the same order
+    binary_maps = np.array([create_binary_map_from_class(k) for k in sorted(class_weights.keys())], dtype=np.uint8)
+
+    logging.debug(f"Creating binary maps with shape {binary_maps.shape}")
+    logging.debug(f"The binary maps order is {sorted(class_weights.keys())}")
+
+    output_tif.parent.mkdir(parents=True, exist_ok=True)
+
+    with rasterio.Env():
+        with rasterio.open(
+            output_tif,
+            "w",
+            driver="GTiff",
+            height=binary_maps.shape[1],
+            width=binary_maps.shape[2],
+            count=binary_maps.shape[0],
+            dtype=rasterio.uint8,
+            crs=crs,
+            transform=rasterio.transform.from_origin(
+                x_min - pixel_size / 2, y_max + pixel_size / 2, pixel_size, pixel_size
+            ),
+        ) as out_file:
+            out_file.write(binary_maps.astype(rasterio.uint8))

coclico/mpla0/mpla0_intrinsic.py

@@ -2,115 +2,22 @@
 import json
 import logging
 from pathlib import Path
-from typing import Dict, Tuple
+from typing import Dict
 
-import laspy
-import numpy as np
-import rasterio
-
-from coclico.metrics.commons import split_composed_class
-
-
-def create_2d_binary_map(
-    xs: np.array,
-    ys: np.array,
-    pixel_size: float,
-    x_min: float,
-    y_max: float,
-    nb_pixels: Tuple[int, int] = (1000, 1000),
-) -> np.array:
-    """Create 2d binary occupancy map from points coordinates:
-    boolean 2d map with 1 where at least one point falls in the pixel, 0 everywhere else.
-
-    Args:
-        xs (np.array): vector of x coordinates of all points
-        ys (np.array): vector of y coordinates of all points
-        pixel_size (float): pixel size (in meters) of the output map
-        x_min (float): x coordinate (in meters) of the pixel center of the upper left corner of the map
-        y_max (float): y coordinate (in meters) of the pixel center of the upper left corner of the map
-        nb_pixels (float, optional): number of pixels on each axis in format (x, y). Defaults to (1000, 1000).
-
-    Returns:
-        np.array: Boolean output map
-    """
-    # numpy array is filled with (y, x) instead of (x, y)
-    grid = np.zeros((nb_pixels[1], nb_pixels[0]), dtype=bool)
-
-    for x, y in zip(xs, ys):
-        grid_x = min(int((x - (x_min - pixel_size / 2)) / pixel_size), nb_pixels[0] - 1)  # x_min is left pixel center
-        grid_y = min(int(((y_max + pixel_size / 2) - y) / pixel_size), nb_pixels[1] - 1)  # y_max is upper pixel center
-        grid[grid_y, grid_x] = True
-
-    return grid
+from coclico.metrics.occupancy_map import create_multilayer_2d_occupancy_map
 
 
 def compute_metric_intrinsic(las_file: Path, class_weights: Dict, output_tif: Path, pixel_size: float = 0.5):
-    """Create 2d occupancy map for each classe that is in class_weights keys, and save result in a single output_tif
-    file with one layer per class (the classes are sorted).
+    """Create 2d occupancy map for each class that is in class_weights keys, and save result in a single output_tif
+    file with one layer per class (the classes are sorted alphabetically).
 
     Args:
         las_file (Path): path to the las file on which to generate mpla0 intrinsic metric
         class_weights (Dict): class weights dict (to know for which classes to generate the binary map)
         output_tif (Path): path to output
         pixel_size (float): size of the output raster pixels
     """
-    with laspy.open(las_file) as f:
-        las = f.read()
-        xs, ys = las.x, las.y
-        classifs = las.classification
-        crs = las.header.parse_crs()
-
-    x_min = np.round(np.min(xs) / pixel_size) * pixel_size
-    y_min = np.round(np.min(ys) / pixel_size) * pixel_size
-    x_max = np.round(np.max(xs) / pixel_size) * pixel_size
-    y_max = np.round(np.max(ys) / pixel_size) * pixel_size
-    logging.debug(f"Raster min/max cell centers infered from in las file: ({x_min}, {x_max}), ({y_min}, {y_max})")
-
-    # add 1 to x_pixel and y_pixel because the map should be able to include the extreme points
-    x_pixels = int(np.ceil((x_max - x_min) / pixel_size)) + 1
-    y_pixels = int(np.ceil((y_max - y_min) / pixel_size)) + 1
-    logging.debug(f"Raster number of pixels infered from in las file: ({x_pixels}, {y_pixels})")
-
-    def create_binary_map_from_class(class_key):
-        splitted_class_key = split_composed_class(class_key)
-        splitted_class_key_int = [int(ii) for ii in splitted_class_key]
-
-        map = create_2d_binary_map(
-            xs[np.isin(classifs, splitted_class_key_int)],
-            ys[np.isin(classifs, splitted_class_key_int)],
-            pixel_size,
-            x_min,
-            y_max,
-            (x_pixels, y_pixels),
-        )
-
-        return map
-
-    # get results for classes that are in weights dictionary (merged if necessary) in a 3d array
-    # which represents a raster with 1 layer per class
-    # keys are sorted to make sure that raster layers can be retrieved in the same order
-    binary_maps = np.array([create_binary_map_from_class(k) for k in sorted(class_weights.keys())], dtype=np.uint8)
-
-    logging.debug(f"Creating binary maps with shape {binary_maps.shape}")
-    logging.debug(f"The binary maps order is {sorted(class_weights.keys())}")
-
-    output_tif.parent.mkdir(parents=True, exist_ok=True)
-
-    with rasterio.Env():
-        with rasterio.open(
-            output_tif,
-            "w",
-            driver="GTiff",
-            height=binary_maps.shape[1],
-            width=binary_maps.shape[2],
-            count=binary_maps.shape[0],
-            dtype=rasterio.uint8,
-            crs=crs,
-            transform=rasterio.transform.from_origin(
-                x_min - pixel_size / 2, y_max + pixel_size / 2, pixel_size, pixel_size
-            ),
-        ) as out_file:
-            out_file.write(binary_maps.astype(rasterio.uint8))
+    create_multilayer_2d_occupancy_map(las_file, class_weights, output_tif, pixel_size)
 
 
 def parse_args():

test/metrics/test_occupancy_map.py

@@ -0,0 +1,58 @@
+import logging
+from pathlib import Path
+
+import numpy as np
+import pytest
+import rasterio
+from pdaltools.las_info import las_info_metadata
+
+from coclico.metrics.occupancy_map import create_multilayer_2d_occupancy_map
+
+pytestmark = pytest.mark.docker
+
+TMP_PATH = Path("./tmp/metrics")
+
+
+def test_create_multilayer_2d_occupancy_map(ensure_test1_data):
+    las_file = Path("./data/test1/niv1/tile_splitted_2818_32247.laz")
+    pixel_size = 0.5
+    las_mtd = las_info_metadata(las_file)
+    las_extent = (las_mtd["minx"], las_mtd["miny"], las_mtd["maxx"], las_mtd["maxy"])
+    logging.debug(f"Test create_multilayer_2d_occupancy_map on las with extent {las_extent}")
+    class_weights = dict(
+        {
+            "0": 1,
+            "1": 1,
+            "2": 0,  # simple classes
+            "3_4_5": 1,  # composed class
+            "3 _ 4": 2,  # composed class with spaces
+        }
+    )
+    output_tif = TMP_PATH / "unit_test_occupancy_map.tif"
+    create_multilayer_2d_occupancy_map(las_file, class_weights, output_tif, pixel_size=pixel_size)
+
+    assert output_tif.exists()
+    with rasterio.Env():
+        with rasterio.open(output_tif) as f:
+            output_data = f.read()
+            output_bounds = f.bounds
+
+    # check that las extent is comprised inside tif extent but tif has no border pixels
+    assert (output_bounds.left < las_extent[0]) and (output_bounds.left > las_extent[0] - pixel_size)  # x_min
+    assert (output_bounds.right > las_extent[2]) and (output_bounds.right < las_extent[2] + pixel_size)  # x_max
+    assert (output_bounds.top > las_extent[3]) and (output_bounds.top < las_extent[3] + pixel_size)  # y_max
+    assert (output_bounds.bottom < las_extent[1]) and (output_bounds.bottom > las_extent[1] - pixel_size)  # y_min
+
+    # Check number of pixels
+    x_las_extent = las_extent[2] - las_extent[0]
+    y_las_extent = las_extent[3] - las_extent[1]
+
+    assert output_data.shape[1] == np.ceil(y_las_extent / pixel_size)
+    assert output_data.shape[2] == np.ceil(x_las_extent / pixel_size)
+
+    assert output_data.shape[0] == len(class_weights.keys())
+    # input las does not contain points with class 0 so output layer should contain only zeroes
+    assert np.all(output_data[0, :, :] == 0)
+    # all other classes have data, so their layers should not contain only zeroes
+    for ii in range(1, len(class_weights.keys())):
+        assert np.any(output_data[ii, :, :] == 1)

test/mpla0/test_mpla0_intrinsic.py

@@ -4,9 +4,7 @@
 import subprocess as sp
 from pathlib import Path
 
-import numpy as np
 import pytest
-import rasterio
 from pdaltools.las_info import las_info_metadata
 
 from coclico.mpla0 import mpla0_intrinsic
@@ -40,23 +38,6 @@ def test_compute_metric_intrinsic(ensure_test1_data):
     mpla0_intrinsic.compute_metric_intrinsic(las_file, class_weights, output_tif, pixel_size=pixel_size)
 
     assert output_tif.exists()
-    with rasterio.Env():
-        with rasterio.open(output_tif) as f:
-            output_data = f.read()
-            output_bounds = f.bounds
-
-    # check that las extent is comprised inside tif extent but tif has no border pixels
-    assert (output_bounds.left < las_extent[0]) and (output_bounds.left > las_extent[0] - pixel_size)  # x_min
-    assert (output_bounds.right > las_extent[2]) and (output_bounds.right < las_extent[2] + pixel_size)  # x_max
-    assert (output_bounds.top > las_extent[3]) and (output_bounds.top < las_extent[3] + pixel_size)  # y_max
-    assert (output_bounds.bottom < las_extent[1]) and (output_bounds.bottom > las_extent[1] - pixel_size)  # y_min
-
-    assert output_data.shape[0] == len(class_weights.keys())
-    # input las does not contain points with class 0 so output layer should contain only zeroes
-    assert np.all(output_data[0, :, :] == 0)
-    # all other classes have data, so their layers should not contain only zeroes
-    for ii in range(1, len(class_weights.keys())):
-        assert np.any(output_data[ii, :, :] == 1)
 
 
 def test_run_main(ensure_test1_data):