Merge pull request #70 from IGNF/fix-pcd-origin

Get tile origin from various point cloud data types

CHANGELOG.md

@@ -1,3 +1,6 @@
+- Add tools to get tile origin from various point cloud data types (las file, numpy array, min/max values)
+- Raise more explicit error when looking a tile origin when the data width is smaller than the buffer size
+
 # 1.7.4
 - Color: fix images bbox to prevent in edge cases where points were at the edge of the last pixel
 - Add possibility to remove points of some classes in standardize

pdaltools/las_info.py

@@ -6,6 +6,8 @@
 import osgeo.osr as osr
 import pdal
 
+from pdaltools.pcd_info import infer_tile_origin
+
 osr.UseExceptions()
 
 
@@ -17,13 +19,37 @@ def las_info_metadata(filename: str):
     return metadata
 
 
-def get_bounds_from_header_info(metadata):
+def get_bounds_from_header_info(metadata: Dict) -> Tuple[float, float, float, float]:
+    """Get bounds from metadata that has been extracted previously from the header of a las file
+
+    Args:
+        metadata (str): Dictonary containing metadata from a las file (as extracted with pipeline.quickinfo)
+
+    Returns:
+        Tuple[float, float, float, float]: minx, maxx, miny, maxy
+    """
     bounds = metadata["bounds"]
     minx, maxx, miny, maxy = bounds["minx"], bounds["maxx"], bounds["miny"], bounds["maxy"]
 
     return minx, maxx, miny, maxy
 
 
+def get_tile_origin_using_header_info(filename: str, tile_width: int = 1000) -> Tuple[int, int]:
+    """ "Get las file theoretical origin (xmin, ymax) for a data that originates from a square tesselation/tiling
+    using the tesselation tile width only, directly from its path
+    Args:
+        filename (str): path to the las file
+        tile_width (int, optional): Tesselation tile width (in meters). Defaults to 1000.
+
+    Returns:
+        Tuple[int, int]: (origin_x, origin_y) tile origin coordinates = theoretical (xmin, ymax)
+    """
+    metadata = las_info_metadata(filename)
+    minx, maxx, miny, maxy = get_bounds_from_header_info(metadata)
+
+    return infer_tile_origin(minx, maxx, miny, maxy, tile_width)
+
+
 def get_epsg_from_header_info(metadata):
     if "srs" not in metadata.keys():
         raise RuntimeError("EPSG could not be inferred from metadata: No 'srs' key in metadata.")

pdaltools/pcd_info.py

@@ -5,42 +5,72 @@
 import numpy as np
 
 
+def infer_tile_origin(minx: float, maxx: float, miny: float, maxy: float, tile_width: int) -> Tuple[int, int]:
+    """Get point cloud theoretical origin (xmin, ymax) for a data that originates from a square tesselation/tiling
+    using the tesselation tile width only, based on the min/max values
+
+    Edge values are supposed to be included in the tile
+
+    Args:
+        minx (float): point cloud min x value
+        maxx (float): point cloud max x value
+        miny (float): point cloud min y value
+        maxy (float): point cloud max y value
+        tile_width (int): tile width in meters
+
+    Raises:
+        ValueError: In case the min and max values do not belong to the same tile
+
+    Returns:
+        Tuple[int, int]: (origin_x, origin_y) tile origin coordinates = theoretical (xmin, ymax)
+    """
+
+    minx_tile_index = np.floor(minx / tile_width)
+    maxx_tile_index = np.floor(maxx / tile_width) if maxx % tile_width != 0 else np.floor(maxx / tile_width) - 1
+    miny_tile_index = np.ceil(miny / tile_width) if miny % tile_width != 0 else np.floor(miny / tile_width) + 1
+    maxy_tile_index = np.ceil(maxy / tile_width)
+
+    if maxx_tile_index == minx_tile_index and maxy_tile_index == miny_tile_index:
+        origin_x = minx_tile_index * tile_width
+        origin_y = maxy_tile_index * tile_width
+        return origin_x, origin_y
+    else:
+        raise ValueError(
+            f"Min values (x={minx} and y={miny}) do not belong to the same theoretical tile as"
+            f"max values (x={maxx} and y={maxy})."
+        )
+
+
 def get_pointcloud_origin_from_tile_width(
     points: np.ndarray, tile_width: int = 1000, buffer_size: float = 0
 ) -> Tuple[int, int]:
     """Get point cloud theoretical origin (xmin, ymax) for a data that originates from a square tesselation/tiling
-    using the tesselation tile width only.
+    using the tesselation tile width only, based on the point cloud as a np.ndarray
 
     Edge values are supposed to be included in the tile
 
+    In case buffer_size is provided, the origin will be calculated on an "original" tile, supposing that
+    there has been a buffer added to the input tile.
 
     Args:
         points (np.ndarray): numpy array with the tile points
         tile_width (int, optional): Edge size of the square used for tiling. Defaults to 1000.
         buffer_size (float, optional): Optional buffer around the tile. Defaults to 0.
 
     Raises:
-        ValueError: Raise an error when the bounding box of the tile is not included in a tile
+        ValueError: Raise an error when the initial tile is smaller than the buffer (in this case, we cannot find the
+        origin (it can be either in the buffer or in the tile))
 
     Returns:
         Tuple[int, int]: (origin_x, origin_y) origin coordinates
     """
     # Extract coordinates xmin, xmax, ymin and ymax of the original tile without buffer
-    x_min, y_min = np.min(points[:, :2], axis=0) + buffer_size
-    x_max, y_max = np.max(points[:, :2], axis=0) - buffer_size
-
-    # Calculate the tiles to which x, y bounds belong
-    tile_x_min = np.floor(x_min / tile_width)
-    tile_x_max = np.floor(x_max / tile_width) if x_max % tile_width != 0 else np.floor(x_max / tile_width) - 1
-    tile_y_min = np.ceil(y_min / tile_width) if y_min % tile_width != 0 else np.floor(y_min / tile_width) + 1
-    tile_y_max = np.ceil(y_max / tile_width)
-
-    if not (tile_x_max - tile_x_min) and not (tile_y_max - tile_y_min):
-        origin_x = tile_x_min * tile_width
-        origin_y = tile_y_max * tile_width
-        return origin_x, origin_y
-    else:
+    minx, miny = np.min(points[:, :2], axis=0) + buffer_size
+    maxx, maxy = np.max(points[:, :2], axis=0) - buffer_size
+
+    if maxx < minx or maxy < miny:
         raise ValueError(
-            f"Min values (x={x_min} and y={y_min}) do not belong to the same theoretical tile as"
-            f"max values (x={x_max} and y={y_max})."
+            "Cannot find pointcloud origin as the pointcloud width or height is smaller than buffer width"
         )
+
+    return infer_tile_origin(minx, maxx, miny, maxy, tile_width)

test/test_las_info.py

@@ -40,6 +40,11 @@ def test_get_bounds_from_quickinfo_metadata():
     assert bounds == (INPUT_MINS[0], INPUT_MAXS[0], INPUT_MINS[1], INPUT_MAXS[1])
 
 
+def test_get_tile_origin_using_header_info():
+    origin_x, origin_y = las_info.get_tile_origin_using_header_info(INPUT_FILE, tile_width=TILE_WIDTH)
+    assert (origin_x, origin_y) == (COORD_X * TILE_COORD_SCALE, COORD_Y * TILE_COORD_SCALE)
+
+
 def test_get_epsg_from_quickinfo_metadata_ok():
     metadata = las_info.las_info_metadata(INPUT_FILE)
     assert las_info.get_epsg_from_header_info(metadata) == "2154"

test/test_pcd_info.py

@@ -12,23 +12,41 @@
 
 
 @pytest.mark.parametrize(
-    "input_points, expected_origin",
+    "minx, maxx, miny, maxy, expected_origin",
     [
-        (np.array([[501, 501, 0], [999, 999, 0]]), (0, 1000)),  # points in the second half
-        (np.array([[1, 1, 0], [400, 400, 0]]), (0, 1000)),  # points in the frist half
-        (np.array([[500, 500, 0], [1000, 500, 0]]), (0, 1000)),  # xmax on edge and xmin in the tile
-        (np.array([[0, 500, 0], [20, 500, 0]]), (0, 1000)),  # xmin on edge and xmax in the tile
-        (np.array([[950, 500, 0], [1000, 500, 0]]), (0, 1000)),  # xmax on edge and xmin in the tile
-        (np.array([[500, 980, 0], [500, 1000, 0]]), (0, 1000)),  # ymax on edge and ymin in the tile
-        (np.array([[500, 0, 0], [500, 20, 0]]), (0, 1000)),  # ymin on edge and ymax in the tile
-        (np.array([[0, 0, 0], [1000, 1000, 0]]), (0, 1000)),  # points at each corner
+        (501, 999, 501, 999, (0, 1000)),  # points in the second half
+        (1, 400, 1, 400, (0, 1000)),  # points in the first half
+        (500, 1000, 500, 500, (0, 1000)),  # xmax on edge and xmin in the tile
+        (0, 20, 500, 500, (0, 1000)),  # xmin on edge and xmax in the tile
+        (950, 1000, 500, 500, (0, 1000)),  # xmax on edge and xmin in the tile
+        (500, 500, 980, 1000, (0, 1000)),  # ymax on edge and ymin in the tile
+        (500, 500, 0, 20, (0, 1000)),  # ymin on edge and ymax in the tile
+        (0, 1000, 0, 1000, (0, 1000)),  # points at each corner
     ],
 )
-def test_get_pointcloud_origin_edge_cases(input_points, expected_origin):
-    origin_x, origin_y = pcd_info.get_pointcloud_origin_from_tile_width(points=input_points, tile_width=1000)
+def test_infer_tile_origin_edge_cases(minx, maxx, miny, maxy, expected_origin):
+    origin_x, origin_y = pcd_info.infer_tile_origin(minx, maxx, miny, maxy, tile_width=1000)
     assert (origin_x, origin_y) == expected_origin
 
 
+@pytest.mark.parametrize(
+    "minx, maxx, miny, maxy",
+    [
+        (0, 20, -1, 20),  # ymin slightly outside the tile
+        (-1, 20, 0, 20),  # xmin slightly outside the tile
+        (280, 1000, 980, 1001),  # ymax slightly outside the tile
+        (980, 1001, 980, 1000),  # xmax slightly outside the tile
+        (-1, 1000, 0, 1000),  # xmax on edge but xmin outside the tile
+        (0, 1000, 0, 1001),  # ymin on edge but ymax outside the tile
+        (0, 1001, 0, 1000),  # xmin on edge but xmax outside the tile
+        (0, 1000, -1, 1000),  # ymax on edge but ymin outside the tile
+    ],
+)
+def test_infer_tile_origin_edge_cases_fail(minx, maxx, miny, maxy):
+    with pytest.raises(ValueError):
+        pcd_info.infer_tile_origin(minx, maxx, miny, maxy, tile_width=1000)
+
+
 @pytest.mark.parametrize(
     "input_points",
     [
@@ -59,3 +77,11 @@ def test_get_pointcloud_origin_on_file():
         points=INPUT_POINTS, tile_width=10, buffer_size=20
     )
     assert (origin_x_2, origin_y_2) == (expected_origin[0] + 20, expected_origin[1] - 20)
+
+
+def test_get_pointcloud_origin_fail_on_buffersize():
+    with pytest.raises(ValueError):
+        # Case when buffer size is bigger than the tile extremities (case not handled)
+        points = np.array([[0, 0, 0], [20, 20, 0]])
+        buffer_size = 30
+        pcd_info.get_pointcloud_origin_from_tile_width(points=points, tile_width=1000, buffer_size=buffer_size)