diff --git a/pdaltools/add_points_in_las.py b/pdaltools/add_points_in_las.py index 46573f3..b47419c 100644 --- a/pdaltools/add_points_in_las.py +++ b/pdaltools/add_points_in_las.py @@ -4,7 +4,7 @@ import numpy as np import pdal -from pdaltools.las_info import get_writer_parameters_from_reader_metadata +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): @@ -12,15 +12,24 @@ def extract_points_from_geo(input_geo: str): 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, values_dimensions: {}): +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] @@ -40,6 +49,7 @@ def parse_args(): 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", @@ -89,5 +99,6 @@ def parse_vars(items): 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/test/test_add_points_in_las.py b/test/test_add_points_in_las.py index 90d26b0..1256406 100644 --- a/test/test_add_points_in_las.py +++ b/test/test_add_points_in_las.py @@ -2,83 +2,71 @@ import os import random as rand import tempfile +import math import pdal import geopandas as gpd -from shapely.geometry import Polygon, Point +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 = 770550 -Ymin = 6277552 +Xmin = 770575 +Ymin = 6277575 Zmin = 20 -Size = 50 - -def build_random_point(): - X = Xmin + rand.uniform(0, 1) * Size - Y = Ymin + rand.uniform(0, 1) * Size - Z = Zmin + rand.uniform(0, 1) * 10 - return Point(X, Y, Z) - -def build_random_geom(nb_geom : int, out_geom_file: str): +Size = 20 - geom=[] - - # add some points - for i in range(nb_geom): - geom.append(build_random_point()) - - # add some polygon: - for i in range(nb_geom): - coordinates = [] - for i in range(4+nb_geom): - coordinates.append(build_random_point()) - polygon = Polygon(coordinates) - geom.append(polygon) +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, only_inside): + geom = [pt_geo] series = gpd.GeoSeries(geom, crs="2154") - series.to_file(out_geom_file) - - # return the number of points - return nb_geom + nb_geom*(4+nb_geom+1) - -@pytest.mark.parametrize("execution_number", range(3)) -def test_extract_points_from_geo(execution_number): - - with tempfile.NamedTemporaryFile(suffix="_geom_tmp.geojson") as geom_file: - nb_points_to_extract = build_random_geom(rand.randint(5,20), geom_file.name) - points = add_points_in_las.extract_points_from_geo(geom_file.name) - assert nb_points_to_extract == len(points) - - -def get_nb_points_from_las(input_las: str, dict_conditions = {}): - pipeline = pdal.Pipeline() | pdal.Reader.las(input_las) - pipeline.execute() - if not dict_conditions: - return len(pipeline.arrays[0]) - return len([e for e in pipeline.arrays[0] if all(e[val] == dict_conditions[val] for val in dict_conditions)]) - - -@pytest.mark.parametrize("execution_number", range(3)) -def test_add_points_in_las(execution_number): with tempfile.NamedTemporaryFile(suffix="_geom_tmp.las") as out_las_file: with tempfile.NamedTemporaryFile(suffix="_geom_tmp.geojson") as geom_file: - nb_points_to_extract = build_random_geom(rand.randint(3, 10), geom_file.name) + 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, added_dimensions) - nb_points_ini = get_nb_points_from_las(INPUT_LAS) - nb_points_to_find = nb_points_ini + nb_points_to_extract - nb_points_end = get_nb_points_from_las(out_las_file.name) - nb_points_end_class = get_nb_points_from_las(out_las_file.name, added_dimensions) - assert nb_points_end == nb_points_to_find - assert nb_points_end_class == nb_points_to_extract + add_points_in_las.add_points_in_las(INPUT_LAS, geom_file.name, out_las_file.name, only_inside, 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, 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, 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, True) + assert len(points_las) == 0