Merge pull request #72 from IGNF/dev

version 1.7.5

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

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

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/_version.py

@@ -1,4 +1,4 @@
-__version__ = "1.7.4"
+__version__ = "1.7.5"
 
 
 if __name__ == "__main__":

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

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.")

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

pdaltools/pcd_info.py

@@ -5,42 +5,72 @@
 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
         tile_width (int, optional): Edge size of the square used for tiling. Defaults to 1000.
         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)

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
 )
 
 

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