Merge pull request #59 from IGNF/get-origin-pointcloud

Add tool to get point cloud origin

CHANGELOG.md

@@ -1,3 +1,5 @@
+# dev
+- Add method to get a point cloud origin
 # 1.7.2
 - Add possibility to select extra dimensions to keep in standardization
 

README.md

@@ -42,6 +42,18 @@ filename = ...
 las_infos.las_info_metadata(filename)
 ```
 
+## Point cloud infos
+
+Misc tools to get information on a point cloud (numpy array). Eg. get expected origin of a point cloud based on a square tiling:
+
+```python
+from pdaltools import pcd_infos
+
+points = ...
+pcd_infos.get_pointcloud_origin_from_tile_width(points, tile_width=1000)
+```
+
+
 ## Stitching
 
 * [las_clip.py](pdaltools/las_clip.py): crop a LAS file using 2d bounding box

environment.yml

@@ -16,6 +16,8 @@ dependencies:
   - pytest
   - laspy
   - requests-mock
-  - build
   - twine
   - tqdm
+  - pip
+  - pip:
+    - build # Installed via pip after issues when creating the environment ("build does not exist")

pdaltools/pcd_info.py

@@ -0,0 +1,46 @@
+"""Tools to get information from a point cloud (points as a numpy array)"""
+
+from typing import Tuple
+
+import numpy as np
+
+
+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.
+
+    Edge values are supposed to be included in the 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
+
+    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:
+        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})."
+        )

test/test_pcd_info.py

@@ -0,0 +1,61 @@
+import os
+
+import laspy
+import numpy as np
+import pytest
+
+from pdaltools import pcd_info
+
+TEST_PATH = os.path.dirname(os.path.abspath(__file__))
+TMP_PATH = os.path.join(TEST_PATH, "tmp")
+DATA_PATH = os.path.join(TEST_PATH, "data")
+
+
+@pytest.mark.parametrize(
+    "input_points, 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
+    ],
+)
+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)
+    assert (origin_x, origin_y) == expected_origin
+
+
+@pytest.mark.parametrize(
+    "input_points",
+    [
+        (np.array([[0, -1, 0], [20, 20, 0]])),  # ymin slightly outside the tile
+        (np.array([[-1, 0, 0], [20, 20, 0]])),  # xmin slightly outside the tile
+        (np.array([[980, 980, 0], [1000, 1001, 0]])),  # ymax slightly outside the tile
+        (np.array([[980, 980, 0], [1001, 1000, 0]])),  # xmax slightly outside the tile
+        (np.array([[-1, 0, 0], [1000, 1000, 0]])),  # xmax on edge but xmin outside the tile
+        (np.array([[0, 0, 0], [1000, 1001, 0]])),  # ymin on edge but ymax outside the tile
+        (np.array([[0, 0, 0], [1001, 1000, 0]])),  # xmin on edge but xmax outside the tile
+        (np.array([[0, -1, 0], [1000, 1000, 0]])),  # ymax on edge but ymin outside the tile
+    ],
+)
+def test_get_pointcloud_origin_edge_cases_fail(input_points):
+    with pytest.raises(ValueError):
+        pcd_info.get_pointcloud_origin_from_tile_width(points=input_points, tile_width=1000)
+
+
+def test_get_pointcloud_origin_on_file():
+    input_las = os.path.join(DATA_PATH, "test_data_77055_627760_LA93_IGN69.laz")
+    expected_origin = (770550, 6277600)
+    LAS = laspy.read(input_las)
+    INPUT_POINTS = np.vstack((LAS.x, LAS.y, LAS.z)).transpose()
+
+    origin_x, origin_y = pcd_info.get_pointcloud_origin_from_tile_width(points=INPUT_POINTS, tile_width=50)
+    assert (origin_x, origin_y) == expected_origin
+    origin_x_2, origin_y_2 = pcd_info.get_pointcloud_origin_from_tile_width(
+        points=INPUT_POINTS, tile_width=10, buffer_size=20
+    )
+    assert (origin_x_2, origin_y_2) == (expected_origin[0] + 20, expected_origin[1] - 20)