diff --git a/CHANGELOG.md b/CHANGELOG.md index b9ffd94..3ce0004 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/pdaltools/las_info.py b/pdaltools/las_info.py index 50a7f4c..3038260 100644 --- a/pdaltools/las_info.py +++ b/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.") diff --git a/pdaltools/pcd_info.py b/pdaltools/pcd_info.py index afdb069..3b40b51 100644 --- a/pdaltools/pcd_info.py +++ b/pdaltools/pcd_info.py @@ -5,14 +5,52 @@ 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 @@ -20,27 +58,19 @@ def get_pointcloud_origin_from_tile_width( 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) diff --git a/test/test_las_info.py b/test/test_las_info.py index d93d8fe..4bcaa8a 100644 --- a/test/test_las_info.py +++ b/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" diff --git a/test/test_pcd_info.py b/test/test_pcd_info.py index 2a08f99..f7687af 100644 --- a/test/test_pcd_info.py +++ b/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)