Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

version 1.7.5 #72

Merged
merged 12 commits into from
Nov 27, 2024
5 changes: 5 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
@@ -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
Expand Down
5 changes: 5 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
1 change: 1 addition & 0 deletions environment.yml
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@ dependencies:
- requests
- gdal
- lastools
- geopandas
# --------- dev dep --------- #
- pre-commit # hooks for applying linters on commit
- black # code formatting
Expand Down
2 changes: 1 addition & 1 deletion pdaltools/_version.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
__version__ = "1.7.4"
__version__ = "1.7.5"


if __name__ == "__main__":
Expand Down
104 changes: 104 additions & 0 deletions pdaltools/add_points_in_las.py
Original file line number Diff line number Diff line change
@@ -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,
)
28 changes: 27 additions & 1 deletion pdaltools/las_info.py
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,8 @@
import osgeo.osr as osr
import pdal

from pdaltools.pcd_info import infer_tile_origin

osr.UseExceptions()


Expand All @@ -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.")
Expand Down
1 change: 0 additions & 1 deletion pdaltools/las_remove_dimensions.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,4 @@
import argparse
import os

import pdal
from pdaltools.las_info import get_writer_parameters_from_reader_metadata
Expand Down
66 changes: 48 additions & 18 deletions pdaltools/pcd_info.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
2 changes: 1 addition & 1 deletion pdaltools/standardize_format.py
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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
)


Expand Down
5 changes: 5 additions & 0 deletions script/test/test_add_points_in_las.sh
Original file line number Diff line number Diff line change
@@ -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
10 changes: 10 additions & 0 deletions test/data/add_points/add_points.geojson
Original file line number Diff line number Diff line change
@@ -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. ] ] ] } }
]
}
72 changes: 72 additions & 0 deletions test/test_add_points_in_las.py
Original file line number Diff line number Diff line change
@@ -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
Loading