From e675205d73d0aac300d1399281d17735d1632e41 Mon Sep 17 00:00:00 2001 From: Antoine Lavenant Date: Fri, 18 Oct 2024 18:00:00 +0200 Subject: [PATCH 1/4] add fn for add points to las --- environment.yml | 1 + pdaltools/add_points_in_las.py | 85 +++++++++++++++++++++++++ pdaltools/las_remove_dimensions.py | 1 - script/test/test_add_points_in_las.sh | 5 ++ test/data/add_points/add_points.geojson | 10 +++ test/test_add_points_in_las.py | 84 ++++++++++++++++++++++++ 6 files changed, 185 insertions(+), 1 deletion(-) create mode 100644 pdaltools/add_points_in_las.py create mode 100644 script/test/test_add_points_in_las.sh create mode 100644 test/data/add_points/add_points.geojson create mode 100644 test/test_add_points_in_las.py 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/add_points_in_las.py b/pdaltools/add_points_in_las.py new file mode 100644 index 0000000..fe5a74d --- /dev/null +++ b/pdaltools/add_points_in_las.py @@ -0,0 +1,85 @@ +import argparse + +import geopandas +import numpy as np +import pdal + +from pdaltools.las_info import get_writer_parameters_from_reader_metadata + +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 add_points_in_las(input_las: str, input_geo: str, output_las: str, 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()) + + for i in points_geo.index: + 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("--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, + values_dimensions=added_dimensions, + ) 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/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..90d26b0 --- /dev/null +++ b/test/test_add_points_in_las.py @@ -0,0 +1,84 @@ +import pytest +import os +import random as rand +import tempfile + +import pdal + +import geopandas as gpd +from shapely.geometry import Polygon, Point + +from pdaltools import add_points_in_las + + + +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 +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): + + 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) + + 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) + 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 + + From 462ecb829489b03d473b95020ab2eaa239b82144 Mon Sep 17 00:00:00 2001 From: Antoine Lavenant Date: Fri, 18 Oct 2024 18:00:50 +0200 Subject: [PATCH 2/4] black --- pdaltools/add_points_in_las.py | 40 ++++++++++++++++++++-------------- 1 file changed, 24 insertions(+), 16 deletions(-) diff --git a/pdaltools/add_points_in_las.py b/pdaltools/add_points_in_las.py index fe5a74d..46573f3 100644 --- a/pdaltools/add_points_in_las.py +++ b/pdaltools/add_points_in_las.py @@ -6,13 +6,13 @@ from pdaltools.las_info import get_writer_parameters_from_reader_metadata + 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 add_points_in_las(input_las: str, input_geo: str, output_las: str, values_dimensions: {}): points_geo = extract_points_from_geo(input_geo) pipeline = pdal.Pipeline() | pdal.Reader.las(input_las) @@ -22,9 +22,9 @@ def add_points_in_las(input_las: str, input_geo: str, output_las: str, values_di for i in points_geo.index: 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] + 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) @@ -40,32 +40,39 @@ 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("--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.") + 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 + if value is None: + return False + try: + nature(value) + return True + except: + return False def parse_var(s): - items = s.split('=') + items = s.split("=") key = items[0].strip() if len(items) > 1: - value = '='.join(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: @@ -74,6 +81,7 @@ def parse_vars(items): d[key] = value return d + if __name__ == "__main__": args = parse_args() added_dimensions = parse_vars(args.dimensions) From 3701e29862e50d2e89de7a7b129504e920a5d205 Mon Sep 17 00:00:00 2001 From: Antoine Lavenant Date: Mon, 18 Nov 2024 16:28:24 +0100 Subject: [PATCH 3/4] ajout d'un test sur l'ajout de point hors du las + fix tests --- pdaltools/add_points_in_las.py | 15 ++++- test/test_add_points_in_las.py | 104 +++++++++++++++------------------ 2 files changed, 59 insertions(+), 60 deletions(-) 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 From d06f81d9f8a5bd59911abc949ac93daf9b8b5c6b Mon Sep 17 00:00:00 2001 From: Antoine Lavenant Date: Mon, 18 Nov 2024 16:53:28 +0100 Subject: [PATCH 4/4] maj nom de variables --- test/test_add_points_in_las.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/test_add_points_in_las.py b/test/test_add_points_in_las.py index 1256406..971bef7 100644 --- a/test/test_add_points_in_las.py +++ b/test/test_add_points_in_las.py @@ -28,7 +28,7 @@ def distance3D(pt_geo, pt_las): numeric_precision, ) -def add_point_in_las(pt_geo, only_inside): +def add_point_in_las(pt_geo, inside_las): geom = [pt_geo] series = gpd.GeoSeries(geom, crs="2154") @@ -37,7 +37,7 @@ def add_point_in_las(pt_geo, only_inside): 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, only_inside, added_dimensions) + 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() @@ -50,7 +50,7 @@ def test_add_point_inside_las(): 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) + 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 @@ -59,7 +59,7 @@ def test_add_point_outside_las_no_control(): 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) + 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 @@ -68,5 +68,5 @@ def test_add_point_outside_las_with_control(): 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) + points_las = add_point_in_las(pt_geo=pt_geo, inside_las=True) assert len(points_las) == 0