-
Notifications
You must be signed in to change notification settings - Fork 1
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #68 from IGNF/add_points_las
Add points in las
- Loading branch information
Showing
6 changed files
with
192 additions
and
1 deletion.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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, | ||
) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 | ||
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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. ] ] ] } } | ||
] | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 |