Add points las

environment.yml

@@ -8,6 +8,7 @@ dependencies:
   - requests
   - gdal
   - lastools
+  - geopandas
 # --------- dev dep --------- #
   - pre-commit # hooks for applying linters on commit
   - black # code formatting

pdaltools/add_points_in_las.py

@@ -0,0 +1,93 @@
+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,
+    )

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

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

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. ] ] ] } }
+]
+}

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
+
+