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

Z on Virtual points #8

Merged
merged 64 commits into from
Jul 18, 2024
Merged
Show file tree
Hide file tree
Changes from 50 commits
Commits
Show all changes
64 commits
Select commit Hold shift + click to select a range
3341e64
add function : crop lidar pointcloud
mdupaysign May 30, 2024
0bc5b80
add test for function : cropping las
mdupaysign May 30, 2024
2b4fa04
add function who calculates quartile and median from np.array
mdupaysign May 30, 2024
4e0ae4f
add test from calculate_stats.py
mdupaysign May 30, 2024
4888491
delete paramaters geom
mdupaysign May 30, 2024
4f87a7f
add function knn from np.array
mdupaysign May 30, 2024
faeecc4
add test for knn
mdupaysign May 30, 2024
3d760fe
modify classe : keep only ground
mdupaysign May 31, 2024
2865602
refacto function calculate_stat.py
mdupaysign Jun 3, 2024
6860779
refacto function knn and test
mdupaysign Jun 3, 2024
b90d069
create function and test: calculate Z value along skeleton
mdupaysign Jun 3, 2024
365db60
rectify test: (output == 0).all() == false
mdupaysign Jun 3, 2024
52c5541
rectify test knn distance
mdupaysign Jun 3, 2024
ab9ad1c
refacto configs with function create virtual point
mdupaysign Jun 6, 2024
cc5c2cd
refacto function crop_las
mdupaysign Jun 6, 2024
dbadb2d
refacto script for lauching steps
mdupaysign Jun 6, 2024
8819752
refacto script las_around_point for lauching
mdupaysign Jun 6, 2024
06f8e6b
refacto test for las around point
mdupaysign Jun 6, 2024
0be63c9
add new function create mask hydro buffer
mdupaysign Jun 6, 2024
b01963d
add function for lauching create virtual point without rectify Z
mdupaysign Jun 6, 2024
12652d2
add test for point along skeleton
mdupaysign Jun 6, 2024
f8d0880
add test for lauching virtual point without rectify Z
mdupaysign Jun 6, 2024
b79e75e
add test for creating mask hydro buffer
mdupaysign Jun 6, 2024
d65f4c9
add function identify upstream/downstream
mdupaysign Jun 7, 2024
eb89910
add function sorted geospatial points by skeleton
mdupaysign Jun 7, 2024
a78d9d3
refacto config.yaml with under-categories
mdupaysign Jun 13, 2024
47e9d25
refacto gitignore: delete all logs
mdupaysign Jun 13, 2024
592eefc
refacto function step 1: import mask hydro and apply buffer
mdupaysign Jun 13, 2024
98c0f21
delete fonction apply_buffer: not only use
mdupaysign Jun 13, 2024
915c4e4
delete fucntion sort_geospatial_points
mdupaysign Jun 13, 2024
82402f2
rename function main_one_tile to extract_points_around_skeleton_point…
mdupaysign Jun 13, 2024
574ca49
add function for extract bounding box with pdal
mdupaysign Jun 14, 2024
2f922d4
add test for utils pdal
mdupaysign Jun 14, 2024
930d731
refacto fucntion las_around_point: N meter default
mdupaysign Jun 14, 2024
d1a0cfa
add fucntion: clip points with bounding box
mdupaysign Jun 14, 2024
6cec9c8
add function: clip point with bounding box
mdupaysign Jun 14, 2024
0b4e98d
refacto function: main create virtual point
mdupaysign Jun 14, 2024
1abcbaf
update
mdupaysign Jul 10, 2024
ac94d3c
refacto environnement
mdupaysign Jul 10, 2024
cf9b8f8
refacto environnement with python 3.11.8
mdupaysign Jul 10, 2024
a44a1d7
refacto environnement with python 3.11
mdupaysign Jul 10, 2024
ea985b0
force geopandas < 1.*
mdupaysign Jul 11, 2024
3fb6f7d
force geopandas==0.*
mdupaysign Jul 11, 2024
5fcc443
rename function crop_las.py, then refacto other related functions
mdupaysign Jul 11, 2024
55e6524
delete function utils_pdal.py for using ign-pdal-tools
mdupaysign Jul 11, 2024
2c0e6a2
rectify function's comments
mdupaysign Jul 11, 2024
8b57015
rectify function's comments v2
mdupaysign Jul 11, 2024
fca39eb
rectify function's comments v3
mdupaysign Jul 11, 2024
927152c
research KNN from 2D points_skeleton NOT 3D
mdupaysign Jul 11, 2024
d4cc007
refacto severals scripts for lauching create virtual point: step3 and 4
mdupaysign Jul 11, 2024
b2091aa
modify filter.range for script crop_las.py and configs (keep_classes_…
mdupaysign Jul 15, 2024
2b3f6ef
create folder (scripts) for putting severals examples scripts
mdupaysign Jul 15, 2024
6b6dec0
modify test (crop_las.py) to List[int]
mdupaysign Jul 15, 2024
ac84533
add test (check dimension) from test_knn_distance.py
mdupaysign Jul 15, 2024
1ccd79d
rename function calculate_quartile to calculate_percentile
mdupaysign Jul 15, 2024
d34e15e
add test (check dimension) from test_knn_distance.py : V2
mdupaysign Jul 15, 2024
fa4befa
add test (check dimension) from test_knn_distance.py : V2bis
mdupaysign Jul 15, 2024
acf7ae7
add test (check dimension) from test_knn_distance.py : V3
mdupaysign Jul 15, 2024
54c35f4
refacto create mask buffer hydro (limit bank + n meters)
mdupaysign Jul 17, 2024
d1b74a0
modify names keep_classes_ground to keep_neighbors_classes
mdupaysign Jul 17, 2024
788b6cb
modify definition : classes
mdupaysign Jul 17, 2024
fe97c9a
modify definition q
mdupaysign Jul 17, 2024
c920aa5
refacto test with classes=[1, 2]
mdupaysign Jul 17, 2024
a0192f6
add test (check dimension) from test_knn_distance.py : term
mdupaysign Jul 17, 2024
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion .gitignore
Original file line number Diff line number Diff line change
@@ -1,3 +1,3 @@
**/__pycache__
tmp
main.log
**.log
45 changes: 28 additions & 17 deletions configs/configs_lidro.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,8 @@ work_dir: ${hydra:runtime.cwd}

io:
input_filename: null
input_mask_hydro: null
input_skeleton: null
input_dir: null
output_dir: null
srid: 2154
Expand All @@ -17,24 +19,33 @@ io:
no_data_value: -9999
tile_size: 1000

raster:
# size for dilatation
dilation_size: 3
mask_generation:
raster:
# size for dilatation
dilation_size: 3
filter:
# Classes to be considered as "non-water"
keep_classes: [0, 1, 2, 3, 4, 5, 6, 17, 64, 65, 66, 67] # All classes
vector:
# Filter water's area (m²)
min_water_area: 150
# Parameters for buffer
buffer_positive: 1
buffer_negative: -1.5 # negative buffer should be bigger than positive buffer to prevent protruding over the banks
# Tolerance from Douglas-Peucker
tolerance: 1

vector:
# Filter water's area (m²)
min_water_area: 150
# Parameters for buffer
buffer_positive: 1
buffer_negative: -1.5 # negative buffer should be bigger than positive buffer to prevent protruding over the banks
# Tolerance from Douglas-Peucker
tolerance: 1
# distance in meters between 2 consecutive points
distance_meters: 2

filter:
# Classes to be considered as "non-water"
keep_classes: [0, 1, 2, 3, 4, 5, 6, 17, 64, 65, 66, 67] # All classes
virtual_point:
filter:
# Keep ground pointcloud between Hydro Mask and Hydro Mask buffer
keep_classes_ground: "[2:2]"
leavauchier marked this conversation as resolved.
Show resolved Hide resolved
vector:
# distance in meters between 2 consecutive points
distance_meters: 2
# buffer for searching the points classification (default. "2") of the input LAS/LAZ file
buffer: 2
# The number of nearest neighbors to find with KNeighbors
k: 20

hydra:
output_subdir: null
Expand Down
4 changes: 3 additions & 1 deletion environment.yml
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@ dependencies:
- scipy
- geojson
- rasterio
- geopandas
- geopandas==0.*
- pyproj
- pdal>=2.6
- python-pdal>=3.2.1
Expand All @@ -24,4 +24,6 @@ dependencies:
- flake8
- isort
- pre-commit
- pip:
- -r requirements.txt

10 changes: 10 additions & 0 deletions example_create_virtual_point_by_tile.sh
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
# For lauching create virtual points
leavauchier marked this conversation as resolved.
Show resolved Hide resolved
python -m lidro.main_create_virtual_point \
io.input_dir=./data/ \
io.input_filename=Semis_2021_0830_6291_LA93_IGN69.laz \
io.input_mask_hydro=./data/merge_mask_hydro/MaskHydro_merge.geosjon \
io.input_skeleton=./data/skeleton_hydro/skeleton_hydro.geojson \
io.output_dir=./tmp/ \



9 changes: 9 additions & 0 deletions example_create_virtual_point_default.sh
Original file line number Diff line number Diff line change
@@ -0,0 +1,9 @@
# Launch hydro mask merging
python -m lidro.main_create_virtual_point \
io.input_dir=./data/ \
io.input_mask_hydro=./data/merge_mask_hydro/MaskHydro_merge.geosjon \
io.input_skeleton=./data/skeleton_hydro/skeleton_hydro.geojson \
io.output_dir=./tmp/ \



2 changes: 1 addition & 1 deletion example_merge_mask_default.sh
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
# For lauching Mask Hydro
# For lauching Mask Hydro merged
leavauchier marked this conversation as resolved.
Show resolved Hide resolved
python -m lidro.main_merge_mask \
io.input_dir=./data/mask_hydro/ \
io.output_dir=./tmp/merge_mask_hydro/ \
Expand Down
35 changes: 35 additions & 0 deletions lidro/create_virtual_point/pointcloud/crop_las.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,35 @@
# -*- coding: utf-8 -*-
""" Filter pointcloud to keep only the ground class """
from typing import List

import numpy as np
import pdal
from shapely.geometry import MultiPolygon


def read_filter_and_crop_pointcloud(input_points: str, geom: MultiPolygon, classes: List[int:int]) -> np.array:
"""Filter pointcloud to keep only the ground class

Args:
input_points (str): Path to the input LAS/LAZ file
geom (MultiPolygon): An array of WKT or GeoJSON 2D MultiPolygon (Mask Hydro with buffer)
classes (list): List of classes to use for the filtering

Returns:
np.array : Numpy array containing point coordinates (X, Y, Z) after filtering and croping
"""
# Crop pointcloud by point
pipeline = (
pdal.Reader.las(filename=input_points, nosrs=True)
| pdal.Filter.range(
limits=f"Classification{classes}",
)
| pdal.Filter.crop(
polygon=geom,
)
)
pipeline.execute()
# extract points
cropped_points = pipeline.arrays[0]

return np.array([cropped_points["X"], cropped_points["Y"], cropped_points["Z"]]).T
37 changes: 37 additions & 0 deletions lidro/create_virtual_point/stats/calculate_stat.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,37 @@
# -*- coding: utf-8 -*-
""" Calculates statistics from np.array (ex. pointcloud) """
import numpy as np


def calculate_quartile(points: np.array, q: int) -> float:
"""
Calculates the quartile's value of Z coordinates

Parameters:
- points (numpy.ndarray): An array of 3D points [x, y, z]
- q (int): Percentage or sequence of percentages for the percentiles to compute.
leavauchier marked this conversation as resolved.
Show resolved Hide resolved
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Comme tu ne testes qu'avec une seule valeur de q dans l'appel de fonction, je pense que tu peux changer ton commentaire pour ne plus parler de séquence de pourcentages (sinon ça implique d'adapter le reste de ton code pour que ça marche aussi) :

Suggested change
- q (int): Percentage or sequence of percentages for the percentiles to compute.
- q (int): Percentage for the percentiles to compute.

Values must be between 0 and 100 inclusive.

Returns:
- float: The quartile of Z coordinates
"""
altitudes = points[:, 2] # Extract the altitude column
n_quartile = np.round(np.percentile(altitudes, q), 2)

return n_quartile


def calculate_median(points: np.array) -> float:
"""
Calculates the median's value of Z coordinates

Parameters:
- points (numpy.ndarray): An array of 3D points [x, y, z]

Returns:
- float: The median of Z coordinates
"""
altitudes = points[:, 2] # Extract the altitude column
n_median = np.round(np.median(altitudes), 2)

return n_median
37 changes: 37 additions & 0 deletions lidro/create_virtual_point/stats/knn_distance.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,37 @@
# -*- coding: utf-8 -*-
""" KNN """
import numpy as np
from sklearn.neighbors import NearestNeighbors


def find_k_nearest_neighbors(points_skeleton: np.array, points_ground_3d: np.array, k: int) -> np.array:
"""Finds the K nearest neighbors of a given point from a list of 3D points

Args:
points_skeleton (np.array): a single point representing the skeleton Hydro
an array of 2D point [x, y]
points_ground_3d (np.array): ground's points from LIDAR.
an array of 3D points [x, y, z]
k (int): the number of nearest neighbors to find

Returns:
numpy.ndarray: an array of the K nearest 3D points [x, y, z]
"""
# Extract 2D coordinates (x, y) from the 3D points
points_ground_2d = points_ground_3d[:, :2]

# Convert point_2d to nump.array
points_skeleton_array = np.array(points_skeleton).reshape(1, -1)

# Initialize and fit the NearestNeighbors model
nbrs = NearestNeighbors(n_neighbors=k, algorithm="auto").fit(points_ground_2d)

# Find the distances and indices of the nearest neighbors
indices = nbrs.kneighbors(points_skeleton_array, return_distance=False)

# Retrieve the points corresponding to the indices
# Use indices to retrieve the closest 3D points
k_nearest_points = points_ground_3d[indices[0]]
# = points_3d[indices.flatten()]

return k_nearest_points
Original file line number Diff line number Diff line change
@@ -0,0 +1,30 @@
# -*- coding: utf-8 -*-
""" Clip Skeleton points by tile (pointcloud)
"""
import geopandas as gpd
from shapely.geometry import box


def clip_points_with_box(gdf_points: gpd.GeoDataFrame, bbox: tuple) -> gpd.GeoDataFrame:
"""Clip skeleton points by tile (bounding box)

Args:
gdf_points (gpd.GeoDataFrame): Points every 2 meters (by default) along skeleton hydro
bbox (tuple): bounding box from tile (pointcloud)

Returns:
gpd.GeoDataframe : Points every 2 meters (by default) along skeleton hydro by tile
"""
# Extract the bounding box limits
xmin = bbox[0][0]
xmax = bbox[0][1]
ymin = bbox[1][0]
ymax = bbox[1][1]

# Create a polygon representing the bounding box
bounding_box = box(xmin, ymin, xmax, ymax)

# Clip points using the bounding box
clipped_points = gdf_points[gdf_points.within(bounding_box)]

return clipped_points
Original file line number Diff line number Diff line change
@@ -0,0 +1,63 @@
# -*- coding: utf-8 -*-
""" Extract points around skeleton by tile
"""
import logging
import os
from typing import List

import geopandas as gpd
from pdaltools.las_info import las_get_xy_bounds
from shapely.geometry import Point

from lidro.create_virtual_point.pointcloud.crop_las import (
read_filter_and_crop_pointcloud,
)
from lidro.create_virtual_point.vectors.clip_points_with_bounding_box import (
clip_points_with_box,
)
from lidro.create_virtual_point.vectors.las_around_point import filter_las_around_point


def extract_points_around_skeleton_points_one_tile(
filename: str,
input_dir: str,
input_mask_hydro_buffer: gpd.GeoDataFrame,
points_skeleton_gdf: gpd.GeoDataFrame,
classes: List[int:int],
k: int,
):
"""Severals steps :
1- Crop filtered pointcloud by Mask Hydro with buffer
2- Extract a Z elevation value along the hydrographic skeleton

Args:
filename (str): filename to the LAS file
input_dir (str): input folder
input_mask_hydro_buffer (gpd.GeoDataFrame): hydro mask with buffer
points_skeleton_gdf (gpd.GeoDataFrame): Points every 2 meters (by default) along skeleton hydro
classes (List): List of classes to use for the filtering
k (int): the number of nearest neighbors to find

Returns:
points_clip (np.array) : Numpy array containing point coordinates (X, Y, Z) after filtering and croping
"""
# Step 1 : Crop filtered pointcloud by Mask Hydro with buffer
input_dir_points = os.path.join(input_dir, "pointcloud")
tilename = os.path.splitext(filename)[0] # filename to the LAS file
input_pointcloud = os.path.join(input_dir_points, filename) # path to the LAS file
logging.info(f"\nCroping filtered pointcloud by Mask Hydro with buffer for tile : {tilename}")
# Croping filtered pointcloud by Mask Hydro with buffer for tile
points_clip = read_filter_and_crop_pointcloud(input_pointcloud, str(input_mask_hydro_buffer), classes)
logging.info(f"\nCropping skeleton points for tile: {tilename}")
# Extract bounding box for clipping points by tile
bbox = las_get_xy_bounds(input_pointcloud)
points_skeleton_clip = clip_points_with_box(points_skeleton_gdf, bbox)
# Create list with Skeleton's points with Z for step 4
points_skeleton_with_z_clip = [
([geom.x, geom.y]) for geom in points_skeleton_clip.geometry if isinstance(geom, Point)
]
# Step 2 : Extract a Z elevation value along the hydrographic skeleton
logging.info(f"\nExtract a Z elevation value along the hydrographic skeleton for tile : {tilename}")
result = filter_las_around_point(points_skeleton_with_z_clip, points_clip, k)

return result
37 changes: 37 additions & 0 deletions lidro/create_virtual_point/vectors/las_around_point.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,37 @@
# -*- coding: utf-8 -*-
""" Extract a Z elevation value every N meters along the hydrographic skeleton
"""
from typing import List

import numpy as np
from shapely.geometry import Point

from lidro.create_virtual_point.stats.calculate_stat import calculate_quartile
from lidro.create_virtual_point.stats.knn_distance import find_k_nearest_neighbors


def filter_las_around_point(points_skeleton: List, points_clip: np.array, k: int) -> List:
"""Extract a Z elevation value every N meters along the hydrographic skeleton

Args:
points_skeleton (list) : points every N meters (by default: 2) along skeleton Hydro
points_clip (np.array): Numpy array containing point coordinates (X, Y, Z) after filtering and croping
k (int): The number of nearest neighbors to find

Returns:
List : Result {'geometry': Point 3D, 'z_q1': points KNN}
"""
# Finds the K nearest neighbors of a given point from a list of 3D points
points_knn_list = [
({"geometry": Point(point), "points_knn": find_k_nearest_neighbors(point, points_clip, k)})
for point in points_skeleton
]

# Calcule Z "Q1" for each points every N meters along skeleton hydro
results = [
({"geometry": p["geometry"], "z_q1": calculate_quartile(p["points_knn"], 10)})
for p in points_knn_list
if not np.all(np.isnan(p["points_knn"]))
]

return results
25 changes: 25 additions & 0 deletions lidro/create_virtual_point/vectors/mask_hydro_with_buffer.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,25 @@
# -*- coding: utf-8 -*-
""" Extract a Z elevation value every N meters along the hydrographic skeleton
"""
import geopandas as gpd
from shapely.geometry import CAP_STYLE


def import_mask_hydro_with_buffer(file_mask: str, buffer: float, crs: str | int) -> gpd.GeoDataFrame:
"""Apply buffer (2 meters by default) from Mask Hydro

Args:
file_mask (str): Path from Mask Hydro
buffer (float): buffer to apply to the input geometry
crs (str | int): Make a CRS from an EPSG code : CRS WKT string

Returns:
gpd.GeoDataFrame : geometry columns are encoded to WKT
"""
# Read Mask Hydro merged
gdf = gpd.read_file(file_mask, crs=crs).unary_union

# Apply buffer (2 meters by default) from Mask Hydro
gdf_buffer = gdf.buffer(buffer, cap_style=CAP_STYLE.square)

return gdf_buffer
4 changes: 2 additions & 2 deletions lidro/main_create_mask.py
Original file line number Diff line number Diff line change
Expand Up @@ -45,8 +45,8 @@ def main(config: DictConfig):
pixel_size = config.io.pixel_size
tile_size = config.io.tile_size
crs = CRS.from_user_input(config.io.srid)
classe = config.filter.keep_classes
dilation_size = config.raster.dilation_size
classe = config.mask_generation.filter.keep_classes
dilation_size = config.mask_generation.raster.dilation_size

def main_on_one_tile(filename):
"""Lauch main.py on one tile
Expand Down
Loading