diff --git a/CHANGELOG.md b/CHANGELOG.md index b9ffd94..6a95ae7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,8 @@ +# 1.7.5 +- 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 +- Add method to add points from vector files (ex : shp, geojson, ...) inside las + # 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/README.md b/README.md index 31dd0cc..e44a1d6 100644 --- a/README.md +++ b/README.md @@ -79,6 +79,11 @@ By default, `xcoord` and `ycoord` are given in kilometers and the shape of the t `readers.las: Global encoding WKT flag not set for point format 6 - 10.` which is due to TerraSolid malformed LAS output for LAS1.4 files with point format 6 to 10. +## Add points in Las + +[add_points_in_las.py](pdaltools/add_points_in_las.py): add points from some vector files (ex: shp, geojson, ...) inside Las. New points will have X,Y and Z coordinates. Other attributes values given by the initial las file are null (ex: classification at 0). These others attributes could be forced by using the '--dimensions/-d' option in the command line (ex : 'add_points_in_las.py -i myLas.las -g myPoints.json -d classification=64' - points will have their classification set to 64). The dimension should be present in the initial las ; this is not allowed to add new dimension. + + # Dev / Build ## Contribute diff --git a/environment.yml b/environment.yml index 7d8ffed..8f1c8c4 100644 --- a/environment.yml +++ b/environment.yml @@ -8,6 +8,7 @@ dependencies: - requests - gdal - lastools + - geopandas # --------- dev dep --------- # - pre-commit # hooks for applying linters on commit - black # code formatting diff --git a/pdaltools/_version.py b/pdaltools/_version.py index 31890a2..d6d2f7b 100644 --- a/pdaltools/_version.py +++ b/pdaltools/_version.py @@ -1,4 +1,4 @@ -__version__ = "1.7.4" +__version__ = "1.7.5" if __name__ == "__main__": diff --git a/pdaltools/add_points_in_las.py b/pdaltools/add_points_in_las.py new file mode 100644 index 0000000..b47419c --- /dev/null +++ b/pdaltools/add_points_in_las.py @@ -0,0 +1,104 @@ +import argparse + +import geopandas +import numpy as np +import pdal + +from pdaltools.las_info import get_writer_parameters_from_reader_metadata, las_info_metadata, get_bounds_from_header_info + + +def extract_points_from_geo(input_geo: str): + file = open(input_geo) + df = geopandas.read_file(file) + return df.get_coordinates(ignore_index=True, include_z=True) + +def point_in_bound(bound_minx, bound_maxx, bound_miny, bound_maxy, pt_x, pt_y): + return pt_x >= bound_minx and pt_x <= bound_maxx and pt_y >= bound_miny and pt_y <= bound_maxy + +def add_points_in_las(input_las: str, input_geo: str, output_las: str, inside_las: bool, values_dimensions: {}): + points_geo = extract_points_from_geo(input_geo) + pipeline = pdal.Pipeline() | pdal.Reader.las(input_las) + pipeline.execute() + points_las = pipeline.arrays[0] + dimensions = list(points_las.dtype.fields.keys()) + + if inside_las: + mtd = las_info_metadata(input_las) + bound_minx, bound_maxx, bound_miny, bound_maxy = get_bounds_from_header_info(mtd) + + for i in points_geo.index: + if inside_las : + if not point_in_bound(bound_minx, bound_maxx, bound_miny, bound_maxy, points_geo["x"][i], points_geo["y"][i]): + continue + pt_las = np.empty(1, dtype=points_las.dtype) + pt_las[0][dimensions.index("X")] = points_geo["x"][i] + pt_las[0][dimensions.index("Y")] = points_geo["y"][i] + pt_las[0][dimensions.index("Z")] = points_geo["z"][i] + for val in values_dimensions: + pt_las[0][dimensions.index(val)] = values_dimensions[val] + points_las = np.append(points_las, pt_las, axis=0) + + params = get_writer_parameters_from_reader_metadata(pipeline.metadata) + pipeline_end = pdal.Pipeline(arrays=[points_las]) + pipeline_end |= pdal.Writer.las(output_las, forward="all", **params) + pipeline_end.execute() + + +def parse_args(): + parser = argparse.ArgumentParser("Add points from geometry file in a las/laz file.") + parser.add_argument("--input_file", "-i", type=str, help="Las/Laz input file") + parser.add_argument("--output_file", "-o", type=str, help="Las/Laz output file.") + parser.add_argument("--input_geo_file", "-g", type=str, help="Geometry input file.") + parser.add_argument("--inside_las", "-l", type=str, help="Keep points only inside the las boundary.") + parser.add_argument( + "--dimensions", + "-d", + metavar="KEY=VALUE", + nargs="+", + help="Set a number of key-value pairs corresponding to value " + "needed in points added in the output las; key should be included in the input las.", + ) + return parser.parse_args() + + +def is_nature(value, nature): + if value is None: + return False + try: + nature(value) + return True + except: + return False + + +def parse_var(s): + items = s.split("=") + key = items[0].strip() + if len(items) > 1: + value = "=".join(items[1:]) + if is_nature(value, int): + value = int(value) + elif is_nature(value, float): + value = float(value) + return (key, value) + + +def parse_vars(items): + d = {} + if items: + for item in items: + key, value = parse_var(item) + d[key] = value + return d + + +if __name__ == "__main__": + args = parse_args() + added_dimensions = parse_vars(args.dimensions) + add_points_in_las( + input_las=args.input_file, + input_geo=args.input_geo_file, + output_las=args.input_file if args.output_file is None else args.output_file, + inside_las=args.inside_las, + values_dimensions=added_dimensions, + ) 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/las_remove_dimensions.py b/pdaltools/las_remove_dimensions.py index ea9e3c3..03931ff 100644 --- a/pdaltools/las_remove_dimensions.py +++ b/pdaltools/las_remove_dimensions.py @@ -1,5 +1,4 @@ import argparse -import os import pdal from pdaltools.las_info import get_writer_parameters_from_reader_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/pdaltools/standardize_format.py b/pdaltools/standardize_format.py index 4b3d71a..e99a794 100644 --- a/pdaltools/standardize_format.py +++ b/pdaltools/standardize_format.py @@ -19,6 +19,7 @@ from pdaltools.unlock_file import copy_and_hack_decorator +# Standard parameters to pass to the pdal writer STANDARD_PARAMETERS = dict( major_version="1", minor_version="4", # Laz format version (pdal always write in 1.x format) @@ -33,7 +34,6 @@ offset_z=0, dataformat_id=6, # No color by default a_srs="EPSG:2154", - class_points_removed=[], # remove points from class ) diff --git a/script/test/test_add_points_in_las.sh b/script/test/test_add_points_in_las.sh new file mode 100644 index 0000000..02dd90e --- /dev/null +++ b/script/test/test_add_points_in_las.sh @@ -0,0 +1,5 @@ +python -u -m pdaltools.add_points_in_las \ + --input_file test/data/test_data_77055_627760_LA93_IGN69.LAZ \ + --output_file test/tmp/addded_cmdline.laz \ + --input_geo_file test/data/add_points/add_points.geojson \ + --dimensions "Classification"=64 "Intensity"=1.1 \ No newline at end of file diff --git a/test/data/add_points/add_points.geojson b/test/data/add_points/add_points.geojson new file mode 100644 index 0000000..129bf7f --- /dev/null +++ b/test/data/add_points/add_points.geojson @@ -0,0 +1,10 @@ +{ +"type": "FeatureCollection", +"name": "add_points", +"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::2154" } }, +"features": [ +{ "type": "Feature", "properties": { "id": 1 }, "geometry": { "type": "MultiLineString", "coordinates": [ [ [ 770558.858057078556158, 6277569.086137800477445, 0. ], [ 770561.225045961793512, 6277574.660014848224819, 0. ], [ 770566.340796128846705, 6277579.241283654235303, 0. ], [ 770566.340796128846705, 6277585.425996542908251, 0. ], [ 770573.747180699137971, 6277589.167366067878902, 0. ], [ 770582.146173510700464, 6277581.455563577823341, 0. ], [ 770587.032860237522982, 6277586.724022705107927, 0. ], [ 770591.308711123419926, 6277590.541746710427105, 0. ], [ 770591.308711123419926, 6277590.541746710427105, 0. ] ] ] } }, +{ "type": "Feature", "properties": { "id": 2 }, "geometry": { "type": "MultiLineString", "coordinates": [ [ [ 770568.402367091737688, 6277561.832462190650403, 0. ], [ 770576.114169582375325, 6277564.123096593655646, 0. ] ] ] } }, +{ "type": "Feature", "properties": { "id": 3 }, "geometry": { "type": "MultiLineString", "coordinates": [ [ [ 770587.18556919763796, 6277575.423559648916125, 0. ], [ 770590.926938722841442, 6277581.91369045805186, 0. ], [ 770593.141218645963818, 6277588.938302627764642, 0. ] ] ] } } +] +} diff --git a/test/test_add_points_in_las.py b/test/test_add_points_in_las.py new file mode 100644 index 0000000..971bef7 --- /dev/null +++ b/test/test_add_points_in_las.py @@ -0,0 +1,72 @@ +import pytest +import os +import random as rand +import tempfile +import math + +import pdal + +import geopandas as gpd +from shapely.geometry import Point + +from pdaltools import add_points_in_las + +numeric_precision = 4 + +TEST_PATH = os.path.dirname(os.path.abspath(__file__)) +INPUT_DIR = os.path.join(TEST_PATH, "data") +INPUT_LAS = os.path.join(INPUT_DIR, "test_data_77055_627760_LA93_IGN69.laz") + +Xmin = 770575 +Ymin = 6277575 +Zmin = 20 +Size = 20 + +def distance3D(pt_geo, pt_las): + return round( + math.sqrt((pt_geo.x - pt_las['X']) ** 2 + (pt_geo.y - pt_las['Y']) ** 2 + (pt_geo.z - pt_las['Z']) ** 2), + numeric_precision, + ) + +def add_point_in_las(pt_geo, inside_las): + geom = [pt_geo] + series = gpd.GeoSeries(geom, crs="2154") + + with tempfile.NamedTemporaryFile(suffix="_geom_tmp.las") as out_las_file: + with tempfile.NamedTemporaryFile(suffix="_geom_tmp.geojson") as geom_file: + series.to_file(geom_file.name) + + added_dimensions = {"Classification":64, "Intensity":1.} + add_points_in_las.add_points_in_las(INPUT_LAS, geom_file.name, out_las_file.name, inside_las, added_dimensions) + + pipeline = pdal.Pipeline() | pdal.Reader.las(out_las_file.name) + pipeline.execute() + points_las = pipeline.arrays[0] + points_las = [e for e in points_las if all(e[val] == added_dimensions[val] for val in added_dimensions)] + return points_las + +def test_add_point_inside_las(): + X = Xmin + rand.uniform(0, 1) * Size + Y = Ymin + rand.uniform(0, 1) * Size + Z = Zmin + rand.uniform(0, 1) * 10 + pt_geo = Point(X, Y, Z) + points_las = add_point_in_las(pt_geo=pt_geo, inside_las=True) + assert len(points_las) == 1 + assert distance3D(pt_geo, points_las[0]) < 1 / numeric_precision + +def test_add_point_outside_las_no_control(): + X = Xmin + rand.uniform(2, 3) * Size + Y = Ymin + rand.uniform(0, 1) * Size + Z = Zmin + rand.uniform(0, 1) * 10 + pt_geo = Point(X, Y, Z) + points_las = add_point_in_las(pt_geo=pt_geo, inside_las=False) + assert len(points_las) == 1 + assert distance3D(pt_geo, points_las[0]) < 1 / numeric_precision + +def test_add_point_outside_las_with_control(): + X = Xmin + rand.uniform(2, 3) * Size + Y = Ymin + rand.uniform(2, 3) * Size + Z = Zmin + rand.uniform(0, 1) * 10 + pt_geo = Point(X, Y, Z) + points_las = add_point_in_las(pt_geo=pt_geo, inside_las=True) + assert len(points_las) == 0 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)