diff --git a/.gitignore b/.gitignore index f711668..b650ae7 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,6 @@ +# QGis generated files +*.copc.laz + # Byte-compiled / optimized / DLL files __pycache__/ *.py[cod] diff --git a/README.md b/README.md index 11db08f..dd9fde1 100644 --- a/README.md +++ b/README.md @@ -45,16 +45,3 @@ DONOR_CLASS_LIST : Défaut [2, 22]. La liste des classes des points du fichier d RECIPIENT_CLASS_LIST : Défaut [2, 3, 9, 17]. La liste des classes des points du fichier receveur qui, s'ils sont absents dans une cellule, justifirons de prendre les points du fichier donneur de la même cellule TILE_SIZE : Défaut 1000. Taille du côté de l'emprise carrée représentée par les fichiers lidar d'entrée PATCH_SIZE : Défaut 1. taille en mètre du côté d'une cellule (doit être un diviseur de TILE_SIZE, soit pour 1000 : 0.25, 0.5, 2, 4, 5, 10, 25...) - -Le script de sélection/découpe de fichier lidar peut être lancé via : -``` -python lidar_selecter.py filepath.DONOR_DIRECTORY=[répertoire_fichiers_donneurs] filepath.RECIPIENT_DIRECTORY=[répertoire_fichiers_receveurs] filepath.SHP_NAME=[nom_shapefile] filepath.SHP_DIRECTORY=[répertoire_shapefile] filepath.CSV_NAME=[nom_fichier_csv] filepath.CSV_DIRECTORY=[répertoire_fichier_csv] filepath.OUTPUT_DIRECTORY=[chemin_de_sortie] -``` - -filepath.DONOR_DIRECTORY: Le répertoire contenant les fichiers lidar donneurs -filepath.RECIPIENT_DIRECTORY: Le répertoire contenant les fichiers lidar receveurs -filepath.SHP_NAME: Le nom du shapefile contenant l'emprise du chantier qui délimite les fichiers lidar qui nous intéressent -filepath.SHP_DIRECTORY: Le répertoire du fichier shapefile -filepath.CSV_NAME: Le nom du fichier csv qui lie les différents fichiers donneurs et receveurs -filepath.CSV_DIRECTORY: Le répertoire du fichier csv -filepath.OUTPUT_DIRECTORY: le répertoire recevant les fichiers lidar découpés \ No newline at end of file diff --git a/configs/configs_patchwork.yaml b/configs/configs_patchwork.yaml index b556c60..7c4a359 100644 --- a/configs/configs_patchwork.yaml +++ b/configs/configs_patchwork.yaml @@ -7,9 +7,9 @@ work_dir: ${hydra:runtime.cwd} # disable ouput directory from being created -hydra: - output_subdir: null - run: +hydra: + output_subdir: null + run: dir: . # disable main.log from being created @@ -19,37 +19,34 @@ defaults: - _self_ filepath: - SHP_NAME: null # name of the shapefile for lidar selecter, to determine the lidar file to select + SHP_NAME: null # name of the shapefile used to match tiles to patch SHP_DIRECTORY: null # path to the directory containing the shapefile - DONOR_DIRECTORY: null # directory containing all potential donor lidar files, for lidar selecter - RECIPIENT_DIRECTORY: null # directory containing all potential donor lidar files, for lidar selecter - OUTPUT_DIRECTORY: null # directory containing all potential donor lidar files, for lidar selecter OUTPUT_DIR: null # directory of the file with added points, from patchwork. OUTPUT_NAME: null # name of the file with added points, from patchwork. - + INPUT_INDICES_MAP_DIR: null INPUT_INDICES_MAP_NAME: null OUTPUT_INDICES_MAP_DIR: null # path to the directory for the indices map reflecting the changes to the recipient, from patchwork OUTPUT_INDICES_MAP_NAME: null # name of the indices map reflecting the changes to the recipient, from patchwork - # INPUT_DIRECTORY: null # directory for input (shapefile) - CSV_NAME: null # name of the csv file that log the lidar files to process with patchwork - CSV_DIRECTORY: null # path to the directory that will contain the csv - - DONOR_NAME: null # name of the donor file for patchwork + RECIPIENT_DIRECTORY: null # directory containing the recipient file for patchwork RECIPIENT_NAME: null # name of the recipient file for patchwork - + # The input shapefile should contain a "nuage_mixa" attrubute for each geometry + # "nuage_mixa" contains the path to the folder containing the files related to a specific donor source. + # Laz/las files from this source are usually contained in a subdirectory of "nuage_mixa" + # path to this subdirectory can be configured using "DONOR_SUBDIRECTORY" + DONOR_SUBDIRECTORY: "data" CRS: 2154 DONOR_CLASS_LIST: [2, 22] RECIPIENT_CLASS_LIST: [2, 6, 9, 17] -RECIPIENT_SUFFIX: "_recipient" TILE_SIZE: 1000 +SHP_X_Y_TO_METER_FACTOR: 1000 # multiplication factor to convert shapefile x, y attributes values to meters PATCH_SIZE: 1 # size of a patch of the grid. Must be a divisor of TILE_SIZE, so for 1000: 0.25, 0.5, 2, 4, 5, 10, 25... NEW_COLUMN: null # If not null, contains the name of the new column NEW_COLUMN_SIZE: 8 # must be 8, 16, 32 or 64 diff --git a/exemples/lidar_selecter_example.sh b/exemples/lidar_selecter_example.sh deleted file mode 100644 index 1f0892b..0000000 --- a/exemples/lidar_selecter_example.sh +++ /dev/null @@ -1,18 +0,0 @@ -# for selecting, cutting and dispatching lidar files for patchwork -python lidar_filepath.py \ -filepath.DONOR_DIRECTORY=[path_to_directory_with_donor_files] \ -filepath.RECIPIENT_DIRECTORY=[path_to_directory_with_recipient_files] \ -filepath.SHP_NAME=[shapefile_name] \ -filepath.SHP_DIRECTORY=[path_to_shapefile_file] \ -filepath.CSV_NAME=[csv_file_name] \ -filepath.CSV_DIRECTORY=[path_to_csv_file] \ -filepath.OUTPUT_DIRECTORY=[output_directory_path] - -# filepath.DONOR_DIRECTORY: The directory containing all the lidar files that could provide points -# filepath.RECIPIENT_DIRECTORY: The directory containing all the lidar files that could receive points -# filepath.SHP_NAME: the name of the shapefile defining the area used to select the lidar files -# filepath.SHP_DIRECTORY: the directory of the shapefile -# filepath.CSV_NAME: the name of the csv file tin which we link donor and recipient files -# filepath.CSV_DIRECTORY: the directory of the csv file -# filepath.OUTPUT_DIRECTORY: the directory to put all the cut lidar files - diff --git a/patchwork/constants.py b/patchwork/constants.py index 260b86b..eece324 100644 --- a/patchwork/constants.py +++ b/patchwork/constants.py @@ -2,9 +2,4 @@ CLASSIFICATION_STR = "classification" PATCH_X_STR = "patch_x" PATCH_Y_STR = "patch_y" -DONOR_SUBDIRECTORY_NAME = "donor" -RECIPIENT_SUBDIRECTORY_NAME = "recipient" - -COORDINATES_KEY = "coordinates" -DONOR_FILE_KEY = "donor_file" -RECIPIENT_FILE_KEY = "recipient_file" +RECIPIENT_SUFFIX = "_recipient" diff --git a/patchwork/lidar_selecter.py b/patchwork/lidar_selecter.py deleted file mode 100644 index 192bd88..0000000 --- a/patchwork/lidar_selecter.py +++ /dev/null @@ -1,153 +0,0 @@ -import os -import pathlib -import timeit - -import geopandas as gpd -import hydra -import laspy -import numpy as np -from laspy import ScaleAwarePointRecord -from loguru import logger -from omegaconf import DictConfig -from pandas import DataFrame -from shapely import box -from shapely.geometry import MultiPolygon -from shapely.vectorized import contains - -import patchwork.constants as c -from patchwork.tools import crop_tile, get_tile_origin_from_pointcloud, identify_bounds - - -@hydra.main(config_path="../configs/", config_name="configs_patchwork.yaml", version_base="1.2") -def patchwork_dispatcher(config: DictConfig): - data = {c.COORDINATES_KEY: [], c.DONOR_FILE_KEY: [], c.RECIPIENT_FILE_KEY: []} - df_result = DataFrame(data=data) - # preparing donor files: - select_lidar( - config, - config.filepath.DONOR_DIRECTORY, - config.filepath.OUTPUT_DIRECTORY, - c.DONOR_SUBDIRECTORY_NAME, - df_result, - c.DONOR_FILE_KEY, - True, - ) - # preparing recipient files: - select_lidar( - config, - config.filepath.RECIPIENT_DIRECTORY, - config.filepath.OUTPUT_DIRECTORY, - c.RECIPIENT_SUBDIRECTORY_NAME, - df_result, - c.RECIPIENT_FILE_KEY, - False, - ) - - pathlib.Path(config.filepath.CSV_DIRECTORY).mkdir(exist_ok=True) - df_result.to_csv( - os.path.join(config.filepath.CSV_DIRECTORY, config.filepath.CSV_NAME), index=False, encoding="utf-8" - ) - - -def cut_lidar(las_points: ScaleAwarePointRecord, shapefile_geometry: MultiPolygon) -> ScaleAwarePointRecord: - shapefile_contains_mask = contains(shapefile_geometry, np.array(las_points["x"]), np.array(las_points["y"])) - return las_points[shapefile_contains_mask] - - -def update_df_result(df_result: DataFrame, df_key: str, corner_string: str, file_path: str): - # corner_string not yet in df_result - if corner_string not in list(df_result[c.COORDINATES_KEY]): - new_row = {c.COORDINATES_KEY: corner_string, c.DONOR_FILE_KEY: "", c.RECIPIENT_FILE_KEY: ""} - new_row[df_key] = file_path - df_result.loc[len(df_result)] = new_row - return df_result - - # corner_string already in df_result - df_result.loc[df_result[c.COORDINATES_KEY] == corner_string, df_key] = file_path - return df_result - - -def select_lidar( - config: DictConfig, - input_directory: str, - output_directory: str, - subdirectory_name: str, - df_result: DataFrame, - df_key: str, - to_be_cut: bool, -): - """ - Walk the input directory searching for las files, and pick the ones that intersect with the shapefile. - When a las file is half inside the shapfile, it is cut if "to_be_cut" is true, otherwise it kept whole - If a file is cut, the cut file is put in: output_directory/subdirectory_name - Finally, df_result is updated with the path for each file - """ - - worksite = gpd.GeoDataFrame.from_file(os.path.join(config.filepath.SHP_DIRECTORY, config.filepath.SHP_NAME)) - shapefile_geometry = worksite.dissolve().geometry.item() - - time_old = timeit.default_timer() - time_start = time_old - - directory_path = os.path.join(output_directory, subdirectory_name) - # Create output dir only if asked to cut - # Otherwise the input file is intended to be used directly (no copy to the output directory) - if to_be_cut: - pathlib.Path(directory_path).mkdir(parents=True, exist_ok=True) - - for root, _, file_names in os.walk(input_directory): - - for file_name in file_names: - if not file_name.endswith((".las", ".laz")): - continue - - logger.info(f"Processing : {file_name}") - las_path = os.path.join(root, file_name) - with laspy.open(las_path) as las_file: - raw_las_points = las_file.read().points - min_x, max_x, min_y, max_y = identify_bounds(config.TILE_SIZE, raw_las_points) - intersect_area = shapefile_geometry.intersection(box(min_x, min_y, max_x, max_y)).area - # if intersect area == 0, this tile is fully outside the shapefile - if intersect_area == 0: - - time_new = timeit.default_timer() - delta_time = round(time_new - time_old, 2) - logger.info(f"Processed {file_name} (out) in {delta_time} sec") - time_old = time_new - continue - - las_points = crop_tile(config, raw_las_points) - x_corner, y_corner = get_tile_origin_from_pointcloud(config, las_points) - - corner_string = f"{int(x_corner/1000)}_{int(y_corner/1000)}" - # if intersect area == TILE_SIZE², this tile is fully inside the shapefile - if intersect_area >= config.TILE_SIZE * config.TILE_SIZE or not to_be_cut: - - df_result = update_df_result(df_result, df_key, corner_string, las_path) - time_new = timeit.default_timer() - delta_time = round(time_new - time_old, 2) - logger.info(f"Processed {file_name} (in) in {delta_time} sec") - time_old = time_new - continue - - # else, this tile is partially inside the shapefile, we have to cut it - points_in_geometry = cut_lidar(las_points, shapefile_geometry) - - # save the selected points to as file - cut_las_path = os.path.join(directory_path, file_name) - with laspy.open(cut_las_path, mode="w", header=las_file.header) as writer: - writer.write_points(points_in_geometry) - - df_result = update_df_result(df_result, df_key, corner_string, cut_las_path) - time_new = timeit.default_timer() - delta_time = round(time_new - time_old, 2) - logger.info(f"Processed {file_name} (cut) in {delta_time} sec") - time_old = time_new - - time_end = timeit.default_timer() - delta_time = round(time_end - time_start, 2) - logger.info(f"Finished processing in {delta_time} sec") - - -if __name__ == "__main__": - patchwork_dispatcher() diff --git a/patchwork/patchwork.py b/patchwork/patchwork.py index 5ad6df8..a645469 100644 --- a/patchwork/patchwork.py +++ b/patchwork/patchwork.py @@ -1,8 +1,8 @@ import os -from pathlib import Path from shutil import copy2 from typing import List, Tuple +import geopandas as gpd import laspy import numpy as np import pandas as pd @@ -12,7 +12,7 @@ import patchwork.constants as c from patchwork.indices_map import create_indices_map -from patchwork.tools import crop_tile, get_tile_origin_from_pointcloud +from patchwork.shapefile_data_extraction import get_donor_info_from_shapefile def get_selected_classes_points( @@ -68,55 +68,60 @@ def get_type(new_column_size: int): raise ValueError(f"{new_column_size} is not a correct value for NEW_COLUMN_SIZE") -def get_complementary_points(config: DictConfig) -> pd.DataFrame: - donor_dir, donor_name = get_donor_path(config) - donor_file_path = os.path.join(donor_dir, donor_name) - recipient_file_path = os.path.join(config.filepath.RECIPIENT_DIRECTORY, config.filepath.RECIPIENT_NAME) - - with laspy.open(donor_file_path) as donor_file, laspy.open(recipient_file_path) as recipient_file: - raw_donor_points = donor_file.read().points - donor_points = crop_tile(config, raw_donor_points) - raw_recipient_points = recipient_file.read().points - recipient_points = crop_tile(config, raw_recipient_points) +def get_complementary_points( + df_donor_info: gpd.GeoDataFrame, recipient_file_path: str, tile_origin: Tuple[int, int], config: DictConfig +) -> pd.DataFrame: + with laspy.open(recipient_file_path) as recipient_file: + recipient_points = recipient_file.read().points + + df_recipient_points = get_selected_classes_points( + config, tile_origin, recipient_points, config.RECIPIENT_CLASS_LIST, [] + ) + + # set, for each patch of coordinate (patch_x, patch_y), the number of recipient point + # should have no record for when count == 0, therefore "df_recipient_non_empty_patches" list all + # and only the patches with at least a point + # In other words, the next column should be filled with "False" everywhere + df_recipient_non_empty_patches = ( + df_recipient_points.groupby(by=[c.PATCH_X_STR, c.PATCH_Y_STR]).count().classification == 0 + ) + + dfs_donor_points = [] + + for index, row in df_donor_info.iterrows(): + with laspy.open(row["full_path"]) as donor_file: + raw_donor_points = donor_file.read().points + points_loc_gdf = gpd.GeoDataFrame( + geometry=gpd.points_from_xy(raw_donor_points.x, raw_donor_points.y, raw_donor_points.z, crs=config.CRS) + ) + footprint_gdf = gpd.GeoDataFrame(geometry=[row["geometry"]], crs=config.CRS) + points_in_footprint_gdf = points_loc_gdf.sjoin(footprint_gdf, how="inner", predicate="intersects") + donor_points = raw_donor_points[points_in_footprint_gdf.index.values] - # check if both files are on the same area - tile_origin_donor = get_tile_origin_from_pointcloud(config, donor_points) - tile_origin_recipient = get_tile_origin_from_pointcloud(config, recipient_points) - if tile_origin_donor != tile_origin_recipient: - raise ValueError( - f"{donor_file_path} and \ - {recipient_file_path} are not on the same area" + donor_columns = get_field_from_header(donor_file) + dfs_donor_points.append( + get_selected_classes_points(config, tile_origin, donor_points, config.DONOR_CLASS_LIST, donor_columns) ) - donor_columns = get_field_from_header(donor_file) - df_donor_points = get_selected_classes_points( - config, tile_origin_donor, donor_points, config.DONOR_CLASS_LIST, donor_columns - ) - df_recipient_points = get_selected_classes_points( - config, tile_origin_recipient, recipient_points, config.RECIPIENT_CLASS_LIST, [] - ) + if len(df_donor_info.index): + df_donor_points = pd.concat(dfs_donor_points) - # set, for each patch of coordinate (patch_x, patch_y), the number of recipient point - # should have no record for when count == 0, therefore "df_recipient_non_empty_patches" list all - # and only the patches with at least a point - # In other words, the next column should be filled with "False" everywhere - df_recipient_non_empty_patches = ( - df_recipient_points.groupby(by=[c.PATCH_X_STR, c.PATCH_Y_STR]).count().classification == 0 - ) + else: + df_donor_points = gpd.GeoDataFrame(columns=["x", "y", "z", "patch_x", "patch_y", "classification"]) - # for each (patch_x,patch_y) patch, we join to a donor point the count of recipient points on that patch - # since it's a left join, it keeps all the left record (all the donor points) - # and put a "NaN" if the recipient point count is null (no record) - joined_patches = pd.merge( - df_donor_points, - df_recipient_non_empty_patches, - on=[c.PATCH_X_STR, c.PATCH_Y_STR], - how="left", - suffixes=("", config.RECIPIENT_SUFFIX), - ) + # for each (patch_x,patch_y) patch, we join to a donor point the count of recipient points on that patch + # since it's a left join, it keeps all the left record (all the donor points) + # and put a "NaN" if the recipient point count is null (no record) + joined_patches = pd.merge( + df_donor_points, + df_recipient_non_empty_patches, + on=[c.PATCH_X_STR, c.PATCH_Y_STR], + how="left", + suffixes=("", c.RECIPIENT_SUFFIX), + ) - # only keep donor points in patches where there is no recipient point - return joined_patches.loc[joined_patches[c.CLASSIFICATION_STR + config.RECIPIENT_SUFFIX].isnull()] + # only keep donor points in patches where there is no recipient point + return joined_patches.loc[joined_patches[c.CLASSIFICATION_STR + c.RECIPIENT_SUFFIX].isnull()] def get_field_from_header(las_file: LasReader) -> List[str]: @@ -144,7 +149,7 @@ def append_points(config: DictConfig, extra_points: pd.DataFrame): c.PATCH_X_STR, c.PATCH_Y_STR, c.CLASSIFICATION_STR, - c.CLASSIFICATION_STR + config.RECIPIENT_SUFFIX, + c.CLASSIFICATION_STR + c.RECIPIENT_SUFFIX, ] fields_to_keep = [ @@ -155,9 +160,6 @@ def append_points(config: DictConfig, extra_points: pd.DataFrame): copy2(recipient_filepath, output_filepath) - if len(extra_points) == 0: # if no point to add, the job is done after copying the recipient file - return - # if we want a new column, we start by adding its name if config.NEW_COLUMN: if test_field_exists(recipient_filepath, config.NEW_COLUMN): @@ -167,9 +169,18 @@ def append_points(config: DictConfig, extra_points: pd.DataFrame): ) new_column_type = get_type(config.NEW_COLUMN_SIZE) output_las = laspy.read(output_filepath) - output_las.add_extra_dim(laspy.ExtraBytesParams(name=config.NEW_COLUMN, type=new_column_type)) + output_las.add_extra_dim( + laspy.ExtraBytesParams( + name=config.NEW_COLUMN, + type=new_column_type, + description="Point origin: 0=initial las", + ) + ) output_las.write(output_filepath) + if len(extra_points) == 0: # if no point to add, the job is done after copying the recipient file + return + with laspy.open(output_filepath, mode="a") as output_las: # put in a new table all extra points and their values on the fields we want to keep new_points = laspy.ScaleAwarePointRecord.zeros(extra_points.shape[0], header=output_las.header) @@ -192,51 +203,25 @@ def append_points(config: DictConfig, extra_points: pd.DataFrame): output_las.append_points(new_points) -def get_donor_from_csv(recipient_file_path: str, csv_file_path: str) -> str: - """ - check if there is a donor file, in the csv file, matching the recipient file - return the path to that file if it exists - return "" otherwise - """ - df_csv_data = pd.read_csv(csv_file_path) - donor_file_paths = df_csv_data.loc[df_csv_data[c.RECIPIENT_FILE_KEY] == recipient_file_path, c.DONOR_FILE_KEY] - if len(donor_file_paths) == 1: - return donor_file_paths.iloc[0] - elif len(donor_file_paths) == 0: - return "" - else: - raise RuntimeError( - f"Found more than one donor file associated with recipient file {recipient_file_path}." - "Please check the matching csv file" - ) - - -def get_donor_path(config: DictConfig) -> Tuple[str, str]: - """Return a donor directory and a name: - If there is no csv file provided in config, return DONOR_DIRECTORY and DONOR_NAME - if there is a csv file provided, return DONOR_DIRECTORY and DONOR_NAME matching the given RECIPIENT - if there is a csv file provided but no matching DONOR, return "" twice""" - if config.filepath.CSV_DIRECTORY and config.filepath.CSV_NAME: - csv_file_path = os.path.join(config.filepath.CSV_DIRECTORY, config.filepath.CSV_NAME) - recipient_file_path = os.path.join(config.filepath.RECIPIENT_DIRECTORY, config.filepath.RECIPIENT_NAME) - donor_file_path = get_donor_from_csv(recipient_file_path, csv_file_path) - if not donor_file_path: # if there is no matching donor file, we do nothing - return "", "" - return str(Path(donor_file_path).parent), str(Path(donor_file_path).name) - return config.filepath.DONOR_DIRECTORY, config.filepath.DONOR_NAME - - def patchwork(config: DictConfig): - _, donor_name = get_donor_path(config) recipient_filepath = os.path.join(config.filepath.RECIPIENT_DIRECTORY, config.filepath.RECIPIENT_NAME) - if donor_name: - complementary_bd_points = get_complementary_points(config) - append_points(config, complementary_bd_points) - - else: # if no matching donor, we simply copy the recipient to the output without doing anything - output_filepath = os.path.join(config.filepath.OUTPUT_DIR, config.filepath.OUTPUT_NAME) - copy2(recipient_filepath, output_filepath) - complementary_bd_points = pd.DataFrame() # No points to add + origin_x_meters, origin_y_meters = get_tile_origin_using_header_info(recipient_filepath, config.TILE_SIZE) + x_shapefile = origin_x_meters / config.SHP_X_Y_TO_METER_FACTOR + y_shapefile = origin_y_meters / config.SHP_X_Y_TO_METER_FACTOR + + shapefile_path = os.path.join(config.filepath.SHP_DIRECTORY, config.filepath.SHP_NAME) + donor_info_df = get_donor_info_from_shapefile( + shapefile_path, + x_shapefile, + y_shapefile, + config.filepath.DONOR_SUBDIRECTORY, + ) + + complementary_bd_points = get_complementary_points( + donor_info_df, recipient_filepath, (origin_x_meters, origin_y_meters), config + ) + + append_points(config, complementary_bd_points) corner_x, corner_y = get_tile_origin_using_header_info(filename=recipient_filepath, tile_width=config.TILE_SIZE) create_indices_map(config, complementary_bd_points, corner_x, corner_y) diff --git a/patchwork/shapefile_data_extraction.py b/patchwork/shapefile_data_extraction.py new file mode 100644 index 0000000..02bba63 --- /dev/null +++ b/patchwork/shapefile_data_extraction.py @@ -0,0 +1,78 @@ +import fnmatch +import os + +import geopandas as gpd + + +def get_donor_info_from_shapefile(input_shapefile: str, x: int, y: int, tile_subdirectory: str) -> gpd.GeoDataFrame: + """Retrieve paths to all the donor files associated with a given tile (with origin x, y) from a shapefile. + + The shapefile should contain one geometry per donor file, with attributes: + - x: string for the x coordinate of the tile + (can be expressed in km, the unit is defined using x_y_to_meters_factor) + - y: string for the y coordinate of the tile + (can be expressed in km, the unit is defined using x_y_to_meters_factor) + - nom_coord: string indicating if the coordinates are expected to be found in the filename + - nuage_mixa: path to the directory that contains the donor file + + The filename for each donor is found using pattern matching using {x}_{y} in the {nuage_mixa}/{tile_subdirectory} + directory. + + It is stored in the "full_path" column of the output geodataframe + + Args: + input_shapefile (str): Shapefile describing donor files + x (int): x coordinate of the tile for which to get the donors + (in the same unit as in the shapefile, usually km) + y (int): y coordinate of the tile for which to get the donors + (in the same unit as in the shapefile, usually km) + tile_subdirectory (str): subdirectory of "nuage_mixa" in which the donor files are stored + + Raises: + NotImplementedError: if nom_coord is false (case not handled) + FileNotFoundError: if path {nuage_mixa}/{tile_subdirectory} does not exist or is not a directory + FileNotFoundError: if there is no file corresponding to coordinates {x}_{y} in {nuage_mixa}/{tile_subdirectory} + RuntimeError: if there is several file corresponding to coordinates {x}_{y} in {nuage_mixa}/{tile_subdirectory} + + Returns: + gpd.GeoDataFrame: geodataframe with columns ["x", "y", "full_path", "geometry"] for each donor file for the + x, y tile + """ + gdf = gpd.GeoDataFrame.from_file(input_shapefile, encoding="utf-8") + gdf = gdf[(gdf["x"].astype(int) == x) & (gdf["y"].astype(int) == y)] + + if not gdf["nom_coord"].isin(["oui", "yes", "true"]).all(): + unsupported_geometries = gdf[(~gdf["nom_coord"].isin(["oui", "yes", "true"]))] + raise NotImplementedError( + "Patchwork currently works only if coordinates are contained in the las/laz names, " + "which is an information store in the 'nom_coord' attribute). " + f"The following geometries do not match this requirement: {unsupported_geometries}" + ) + + if len(gdf.index): + + def find_las_path_from_geometry_attributes(x: int, y: int, path_root: str): + tile_directory = os.path.join(path_root, tile_subdirectory) + if not os.path.isdir(tile_directory): + raise FileNotFoundError(f"Directory {tile_directory} not found") + potential_filenames = fnmatch.filter(os.listdir(tile_directory), f"*{x}_{y}*.la[sz]") + if not potential_filenames: + raise FileNotFoundError( + f"Could not match any file with directory {tile_directory} and coords ({x}, {y})" + ) + if len(potential_filenames) > 1: + raise RuntimeError( + f"Found multiple files for directory {tile_directory} and coords ({x}, {y}): {potential_filenames}" + ) + + return os.path.join(tile_directory, potential_filenames[0]) + + gdf["full_path"] = gdf.apply( + lambda row: find_las_path_from_geometry_attributes(row["x"], row["y"], row["nuage_mixa"]), + axis="columns", + ) + + else: + gdf = gpd.GeoDataFrame(columns=["x", "y", "full_path", "geometry"]) + + return gdf[["x", "y", "full_path", "geometry"]] diff --git a/patchwork/tools.py b/patchwork/tools.py deleted file mode 100644 index 313e646..0000000 --- a/patchwork/tools.py +++ /dev/null @@ -1,47 +0,0 @@ -from typing import Tuple - -import numpy as np -from laspy import ScaleAwarePointRecord -from omegaconf import DictConfig - - -def get_tile_origin_from_pointcloud(config: DictConfig, points: np.ndarray | ScaleAwarePointRecord) -> Tuple[int, int]: - """Return the coordinate of the upper left corner of an area defined by the points it contains""" - if not len(points): - raise ValueError("No points to determine the coordinate of the tile") - - x_min = np.min(points["x"], axis=0) - x_max = np.max(points["x"], axis=0) - y_min = np.min(points["y"], axis=0) - y_max = np.max(points["y"], axis=0) - - length_tile_x = x_max / config.TILE_SIZE - np.floor(x_min / config.TILE_SIZE) - length_tile_y = np.ceil(y_max / config.TILE_SIZE) - y_min / config.TILE_SIZE - - if (length_tile_x <= 1) and (length_tile_y <= 1): - origin_x = np.floor(x_min / config.TILE_SIZE) * config.TILE_SIZE - origin_y = np.ceil(y_max / config.TILE_SIZE) * config.TILE_SIZE - return origin_x, origin_y - else: - 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})." - ) - - -def identify_bounds(tile_size: float, points: ScaleAwarePointRecord) -> Tuple[int, int, int, int]: - """Return the bounds of a tile represented by its points""" - gravity_x = np.sum(points["x"]) / len(points) - gravity_y = np.sum(points["y"]) / len(points) - min_x = int(gravity_x / tile_size) * tile_size - max_x = int(gravity_x / tile_size) * tile_size + tile_size - min_y = int(gravity_y / tile_size) * tile_size - max_y = int(gravity_y / tile_size) * tile_size + tile_size - - return min_x, max_x, min_y, max_y - - -def crop_tile(config: DictConfig, points: ScaleAwarePointRecord) -> np.ndarray: - """Crop points to the tile containing the center of gravity""" - min_x, max_x, min_y, max_y = identify_bounds(config.TILE_SIZE, points) - return points[(points["x"] >= min_x) & (points["x"] <= max_x) & (points["y"] >= min_y) & (points["y"] <= max_y)] diff --git a/test/data/aveyron_aval_lidarBD/data/NUALID_1-0_VLIDAVEYRONAVAL_PTS_0673_6362_LAMB93_IGN69_20210319.laz b/test/data/aveyron_aval_lidarBD/data/NUALID_1-0_VLIDAVEYRONAVAL_PTS_0673_6362_LAMB93_IGN69_20210319.laz new file mode 100755 index 0000000..433483b Binary files /dev/null and b/test/data/aveyron_aval_lidarBD/data/NUALID_1-0_VLIDAVEYRONAVAL_PTS_0673_6362_LAMB93_IGN69_20210319.laz differ diff --git a/test/data/aveyron_aval_lidarBD/data/NUALID_1-0_VLIDAVEYRONAVAL_PTS_0673_6363_LAMB93_IGN69_20210319.laz b/test/data/aveyron_aval_lidarBD/data/NUALID_1-0_VLIDAVEYRONAVAL_PTS_0673_6363_LAMB93_IGN69_20210319.laz new file mode 100755 index 0000000..3cb7da6 Binary files /dev/null and b/test/data/aveyron_aval_lidarBD/data/NUALID_1-0_VLIDAVEYRONAVAL_PTS_0673_6363_LAMB93_IGN69_20210319.laz differ diff --git a/test/data/aveyron_lidarBD/data/NUALID_1-0_IAVEY_PTS_0673_6363_LAMB93_IGN69_20170519.laz b/test/data/aveyron_lidarBD/data/NUALID_1-0_IAVEY_PTS_0673_6363_LAMB93_IGN69_20170519.laz new file mode 100755 index 0000000..7fc4efe Binary files /dev/null and b/test/data/aveyron_lidarBD/data/NUALID_1-0_IAVEY_PTS_0673_6363_LAMB93_IGN69_20170519.laz differ diff --git a/test/data/aveyron_lidarBD/data/NUALID_1-0_IAVEY_PTS_0673_6364_LAMB93_IGN69_20170519.laz b/test/data/aveyron_lidarBD/data/NUALID_1-0_IAVEY_PTS_0673_6364_LAMB93_IGN69_20170519.laz new file mode 100755 index 0000000..724c6ea Binary files /dev/null and b/test/data/aveyron_lidarBD/data/NUALID_1-0_IAVEY_PTS_0673_6364_LAMB93_IGN69_20170519.laz differ diff --git a/test/data/donor_infos/donor_info_673_6362_one_donor.csv b/test/data/donor_infos/donor_info_673_6362_one_donor.csv new file mode 100644 index 0000000..46a3c60 --- /dev/null +++ b/test/data/donor_infos/donor_info_673_6362_one_donor.csv @@ -0,0 +1,2 @@ +,x,y,full_path,geometry +0,0673,6362,test/data/aveyron_aval_lidarBD/data/NUALID_1-0_VLIDAVEYRONAVAL_PTS_0673_6362_LAMB93_IGN69_20210319.laz,"POLYGON ((673416.6253332305 6361999.533745996, 673419.0434209384 6361992.529629876, 673422.4205009021 6361983.6106751, 673422.9376225754 6361982.213608164, 673424.9722914887 6361976.588347051, 673428.408799128 6361969.041506745, 673431.3470301083 6361961.863034509, 673434.4470301083 6361953.363034509, 673437.7174257919 6361942.525885733, 673439.7174257919 6361934.225885733, 673442.0387030325 6361920.394091859, 673442.8387030326 6361912.09409186, 673443.2783056729 6361904.582881368, 673443.4783056729 6361894.982881368, 673443.4362854348 6361889.330846949, 673443.1362854347 6361880.930846948, 673441.7460115744 6361867.509308349, 673440.2460115744 6361858.809308349, 673436.8124801891 6361844.919324614, 673434.9122094598 6361839.0664907675, 673431.9397664428 6361828.535549793, 673428.0742113594 6361817.398497735, 673426.5743539362 6361813.781194538, 673424.0400038498 6361807.483718566, 673421.4331362228 6361800.387245581, 673420.0875348861 6361796.393848065, 673419.3621920279 6361793.4579364965, 673417.9158137618 6361786.686256432, 673416.7186339626 6361779.714444661, 673414.2555792663 6361763.528656657, 673412.9410582777 6361753.089813512, 673412.5014155972 6361748.180470245, 673407.1722272978 6361723.742134955, 673404.8722272977 6361717.242134955, 673400.3554326419 6361706.50961203, 673397.555432642 6361700.809612031, 673390.121278595 6361688.126704425, 673386.121278595 6361682.326704424, 673385.8244718217 6361681.898723595, 673380.5244718216 6361674.298723595, 673375.9443055603 6361668.235637708, 673370.1443055603 6361661.135637708, 673365.5957299203 6361655.576220974, 673354.7087760442 6361644.223479534, 673348.3087760443 6361638.523479533, 673334.7998939996 6361628.400169599, 673331.5998939995 6361626.400169599, 673317.3042821205 6361618.993825882, 673302.0053838543 6361613.977636283, 673296.6053838542 6361612.677636283, 673290.3341174779 6361611.378824519, 673283.434117478 6361610.178824519, 673281.3797989732 6361609.843540105, 673275.4797989732 6361608.943540105, 673256.6988339456 6361607.868516623, 673251.2988339455 6361608.068516623, 673245.2900035625 6361608.472536609, 673241.1900035625 6361608.872536609, 673238.2368113322 6361609.205022039, 673233.5368113321 6361609.8050220385, 673211.2965697311 6361615.989709211, 673211.1772233931 6361615.631670197, 673205.477223393 6361617.531670197, 673201.4493518078 6361618.9707150655, 673193.8493518078 6361621.870715065, 673187.685639443 6361624.461906387, 673185.8148332557 6361625.323071139, 673166.7885438359 6361619.813008993, 673163.4885438358 6361619.213008992, 673142.6600948233 6361617.643224554, 673139.2600948233 6361617.743224555, 673121.0667207517 6361619.9585834555, 673117.3667207517 6361620.758583455, 673094.8401913649 6361628.534333716, 673088.0401913648 6361631.834333716, 673084.8497370114 6361633.453789793, 673078.2497370114 6361636.953789793, 673063.1355711425 6361646.8117234465, 673061.2355711425 6361648.3117234465, 673047.7097364344 6361661.21631219, 673042.5381745647 6361668.434605423, 673036.8766997271 6361672.822248423, 673031.2275256416 6361677.10192576, 673024.700569088 6361681.915556218, 673022.0026844952 6361683.808102126, 673019.5452520653 6361684.261781959, 673015.2140493266 6361685.160880431, 673010.0140493267 6361686.360880431, 673004.0911670268 6361687.920188731, 673000 6361689.132386369, 673000 6361959.781600325, 673006.7419231678 6361957.452564301, 673013.7419231678 6361954.152564301, 673025.4895454248 6361947.615298656, 673030.8895454248 6361944.115298656, 673036.2823879762 6361940.362747507, 673042.1823879762 6361935.962747507, 673062.0162194129 6361916.308326753, 673065.8162194128 6361911.308326753, 673074.2022172662 6361898.293260113, 673076.9011235306 6361893.292345565, 673079.0197785955 6361889.629585961, 673085.5481225636 6361888.111602976, 673092.4481225637 6361886.211602976, 673094.6347885539 6361885.582628523, 673100.6347885539 6361883.782628523, 673112.0661495088 6361879.57882088, 673117.7661495088 6361877.07882088, 673129.2493170215 6361871.129130856, 673135.5493170215 6361867.329130856, 673142.023819372 6361863.073347121, 673146.4515360891 6361859.910692323, 673151.9556183502 6361856.213920655, 673159.1198228706 6361850.924487067, 673162.5046997464 6361848.184348644, 673167.2209978788 6361844.919219168, 673170.250069046 6361842.675260993, 673174.0562775683 6361843.005615835, 673183.2562775683 6361843.405615835, 673194.6688863962 6361843.249841329, 673207.3688863963 6361842.349841328, 673209.9309709016 6361842.077831461, 673210.7785278414 6361846.645221637, 673211.6608804307 6361850.8859506715, 673213.4608804308 6361858.685950671, 673215.6717711961 6361866.721868208, 673218.1717711961 6361874.5218682075, 673221.2167560435 6361882.758863936, 673224.9167560434 6361891.558863936, 673227.2432484005 6361896.683529848, 673227.5899128705 6361897.39336662, 673227.4805430147 6361897.599322842, 673225.7448363636 6361901.025093417, 673218.7448363636 6361915.525093417, 673217.3938999902 6361918.442326484, 673210.6938999902 6361933.542326484, 673209.9881472522 6361935.171776519, 673204.4854228874 6361948.192307412, 673197.2869610942 6361961.956399191, 673196.0658923903 6361964.370134191, 673187.1658923903 6361982.570134191, 673183.9948010043 6361989.757205335, 673179.9482648404 6362000, 673416.4561834084 6362000, 673416.6253332305 6361999.533745996))" diff --git a/test/data/donor_infos/donor_info_673_6363_two_donors.csv b/test/data/donor_infos/donor_info_673_6363_two_donors.csv new file mode 100644 index 0000000..ed0a494 --- /dev/null +++ b/test/data/donor_infos/donor_info_673_6363_two_donors.csv @@ -0,0 +1,3 @@ +,x,y,full_path,geometry +1,0673,6363,test/data/aveyron_aval_lidarBD/data/NUALID_1-0_VLIDAVEYRONAVAL_PTS_0673_6363_LAMB93_IGN69_20210319.laz,"POLYGON ((673378.6919298545 6362089.5854969965, 673374.3378244891 6362087.71945184, 673376.8770713893 6362080.9286128385, 673379.3021560904 6362078.001786475, 673384.3398166763 6362071.37926617, 673386.6398166763 6362068.07926617, 673397.6751493448 6362047.46523724, 673399.8751493449 6362041.86523724, 673402.1582665141 6362035.413136791, 673403.2001469469 6362032.113848753, 673404.0895333156 6362029.723622888, 673408.3458935957 6362020.643387624, 673411.2618696775 6362013.765136249, 673415.5618696774 6362002.46513625, 673416.4561834084 6362000, 673179.9482648404 6362000, 673178.528428188 6362003.593961527, 673171.4808721666 6362019.081183401, 673171.2577049187 6362019.575147058, 673165.1577049188 6362033.1751470575, 673163.0061337576 6362038.356668478, 673159.9061337577 6362046.456668478, 673155.7492758363 6362060.203268081, 673153.4492758362 6362070.403268081, 673151.6116265327 6362081.356847392, 673150.7116265327 6362089.456847392, 673150.1337154426 6362103.096526872, 673150.3337154427 6362110.796526873, 673150.869084608 6362118.853312364, 673151.469084608 6362124.4533123635, 673152.8969410876 6362133.684678619, 673154.2969410876 6362140.584678619, 673159.3583686741 6362157.603294792, 673162.058368674 6362164.4032947915, 673165.5572808989 6362172.221359548, 673169.2572808989 6362179.6213595485, 673171.8080805743 6362184.395397145, 673176.3080805743 6362192.295397145, 673177.3637751846 6362194.104410234, 673182.5637751848 6362202.804410234, 673183.3496447372 6362204.09693023, 673187.822760228 6362211.330053151, 673191.2990492805 6362217.0612324, 673198.8491065227 6362227.8391108345, 673203.3491065227 6362233.439110834, 673212.3344827583 6362243.213793103, 673216.5344827583 6362247.213793103, 673224.6568456374 6362254.160636132, 673227.6568456374 6362256.460636131, 673244.8726493459 6362267.081410718, 673251.4726493459 6362270.281410718, 673273.406954215 6362277.918706018, 673278.806954215 6362279.118706018, 673283.6193079123 6362279.972600438, 673284.3711100012 6362280.351770187, 673385.069810581 6362092.462507275, 673380.7481700379 6362090.494235344, 673378.6919298545 6362089.5854969965))" +2,0673,6363,test/data/aveyron_lidarBD/data/NUALID_1-0_IAVEY_PTS_0673_6363_LAMB93_IGN69_20170519.laz,"POLYGON ((673600.5005554566 6362991.191050272, 673604.4005435876 6362985.726480441, 673613.3005435877 6362972.326480441, 673619.3256486757 6362961.954738223, 673627.0256486756 6362946.654738223, 673630.7565951694 6362938.3124308875, 673637.9565951695 6362920.012430888, 673638.8974443333 6362917.524484741, 673649.188045648 6362889.178555664, 673666.0443883506 6362848.149106119, 673688.251700554 6362798.280054153, 673688.4496544013 6362797.8325835485, 673707.6114965341 6362754.229616246, 673720.0517285521 6362727.516869975, 673721.0562962577 6362725.289040453, 673732.3562962576 6362699.389040453, 673732.3987672278 6362699.2915541055, 673740.0987672278 6362681.591554105, 673740.8974813808 6362679.702841186, 673746.1974813808 6362666.802841186, 673747.7799194751 6362662.696441577, 673752.6799194751 6362649.096441577, 673753.8369659113 6362645.694594997, 673760.3369659113 6362625.394594997, 673761.6559912023 6362620.918081461, 673765.071327902 6362608.243387487, 673766.1104851161 6362605.365721355, 673769.2621616236 6362596.86574532, 673771.8196301331 6362588.979896782, 673775.4196301331 6362576.0798967825, 673777.6071038123 6362566.414833677, 673779.2703546465 6362556.897342793, 673780.4211915119 6362550.8775807265, 673781.9179282294 6362539.6056505125, 673782.4623181222 6362532.373041939, 673783.2497592705 6362525.64950598, 673784.2705865578 6362519.451626021, 673784.7488349625 6362516.219543985, 673786.0488349625 6362506.319543985, 673786.848011434 6362490.075870598, 673786.5731307453 6362481.554569251, 673786.386114004 6362470.333564767, 673785.7331022813 6362460.4702649135, 673784.4331022813 6362449.270264913, 673783.7393923833 6362444.36010127, 673781.5393923833 6362431.160101269, 673781.1533481755 6362428.991411329, 673778.6533481755 6362415.79141133, 673777.9818293508 6362412.541670226, 673775.1818293508 6362400.041670226, 673773.8373790372 6362394.727092978, 673771.8052656099 6362387.530024589, 673770.2377437385 6362381.3424382545, 673768.6993716902 6362375.917340326, 673767.3817270729 6362371.724834725, 673767.2052578483 6362364.362171383, 673766.5052578482 6362353.162171383, 673764.534970312 6362338.704140897, 673763.434970312 6362333.504140897, 673755.0427191013 6362309.478640453, 673753.4427191013 6362306.278640453, 673740.4291483532 6362286.512906082, 673737.7291483532 6362283.312906082, 673732.0106781187 6362277.0893218815, 673729.5106781187 6362274.5893218815, 673720.3997794222 6362267.810533993, 673716.236336061 6362262.858426505, 673710.036336061 6362256.958426505, 673702.4940613515 6362250.464778263, 673697.9940613515 6362246.964778263, 673689.1269203793 6362240.806388986, 673681.0269203794 6362235.806388986, 673665.6390676381 6362228.052330913, 673654.6390676381 6362223.652330913, 673646.7714272093 6362220.880004449, 673638.4563186722 6362218.33456306, 673635.226379774 6362217.235009393, 673616.4833804844 6362212.81317721, 673597.2403464939 6362212.0660058325, 673592.040346494 6362212.366005832, 673565.0629515594 6362217.843756304, 673562.3892886421 6362215.931141477, 673560.3668545057 6362214.670783971, 673555.5640679651 6362210.849211886, 673549.7683915592 6362206.569273875, 673541.8275774383 6362201.136085266, 673538.418608994 6362198.609438066, 673531.8199992727 6362190.163530103, 673527.0199992727 6362185.063530102, 673517.0869249435 6362175.848892798, 673510.2869249436 6362170.348892798, 673504.5792661656 6362166.06018332, 673496.7471376127 6362160.601427057, 673491.6397741912 6362156.645018772, 673486.0953029888 6362152.645600809, 673477.5953029888 6362146.945600809, 673475.851754813 6362145.802564453, 673465.7710174425 6362139.343062838, 673460.5895454248 6362135.984701345, 673459.2281429314 6362135.206370842, 673456.3428777345 6362132.926710895, 673450.8428777345 6362129.026710895, 673441.0534651963 6362122.902426016, 673433.7534651962 6362118.902426016, 673431.3316647599 6362117.618253718, 673423.53166476 6362113.618253718, 673415.648086927 6362110.211464346, 673415.2052081259 6362109.91809245, 673411.5577961705 6362107.608064878, 673407.5238193724 6362104.726652879, 673390.8481700378 6362095.094235344, 673385.069810581 6362092.462507275, 673284.3711100012 6362280.351770187, 673292.7683469116 6362284.586898368, 673293.4339211908 6362284.919479194, 673300.38352466 6362288.359876951, 673307.2504188012 6362292.382905842, 673311.2568066395 6362294.608367694, 673319.5593744337 6362298.974373172, 673323.6056887265 6362301.503319605, 673330.4498581651 6362306.1475774385, 673333.6700908738 6362308.243530808, 673342.8669971749 6362313.981050335, 673358.763078958 6362324.807570901, 673374.8089735393 6362336.281730154, 673375.3316429212 6362336.67373219, 673376.8486647375 6362338.884249694, 673385.2280965266 6362349.481493849, 673389.1280965265 6362353.781493849, 673396.1086675965 6362360.753577906, 673400.3086675965 6362364.553577906, 673407.9956972105 6362370.843326697, 673414.4956972105 6362375.643326697, 673419.4394557144 6362379.0692382, 673427.1394557144 6362384.0692382, 673428.6001059996 6362384.999830401, 673434.1817832218 6362388.4883786645, 673438.7483283577 6362391.718374005, 673445.2733062016 6362396.653934169, 673450.1299803813 6362400.105029437, 673458.2299803813 6362405.505029436, 673460.1947918833 6362406.781907557, 673465.0065595275 6362409.8293603975, 673468.7399620436 6362412.343692705, 673474.4963729484 6362415.94262855, 673480.1963729485 6362419.24262855, 673494.834576585 6362426.199752631, 673498.9942471489 6362427.777558707, 673502.3767515703 6362429.1982105635, 673524.4349406085 6362435.601601384, 673531.5349406084 6362436.801601384, 673535.1965963232 6362437.350953061, 673541.2965963233 6362438.150953061, 673559.6069402691 6362438.855177431, 673562.0313400687 6362445.907976849, 673563.9815075943 6362451.575665096, 673564.7113738022 6362454.3870016, 673565.3847753962 6362458.202943966, 673565.9827252551 6362464.834751492, 673566.6071780241 6362474.7479392, 673566.7816757273 6362479.866538492, 673566.7057805138 6362486.9247933375, 673566.7 6362488, 673566.7 6362493.906379602, 673566.3413981188 6362498.388903118, 673565.6256612801 6362504.830534665, 673565.3135255335 6362506.904007839, 673564.5769034621 6362509.710187159, 673564.4573177648 6362510.170200612, 673562.4484010994 6362517.974069197, 673560.6808605859 6362523.700900461, 673557.9143223667 6362532.593344737, 673557.8476412098 6362532.80853124, 673555.3476412098 6362540.90853124, 673555.0742437382 6362541.809364526, 673546.660739834 6362570.008471018, 673538.7116077441 6362592.626620057, 673528.0875480856 6362618.070796549, 673514.7503218431 6362649.252311696, 673500.5590472363 6362678.284134911, 673500.2514259935 6362678.919119642, 673482.1514259935 6362716.619119642, 673479.6717053219 6362722.217125574, 673461.7717053219 6362766.217125574, 673460.777184575 6362768.760728069, 673446.4771845749 6362806.860728068, 673445.3950835699 6362809.891141348, 673435.8388907175 6362838.077083902, 673430.4465923568 6362852.928003978, 673429.1491351919 6362856.289597541, 673412.2927160816 6362869.792858911, 673407.2927160816 6362874.892858911, 673395.7157430674 6362889.100241028, 673391.8157430675 6362894.900241028, 673383.0654988913 6362910.890688196, 673378.537567546 6362921.324616947, 673371.510237072 6362937.358807607, 673370.0467583142 6362940.879046821, 673363.9467583143 6362956.379046821, 673361.8777181053 6362962.149530196, 673357.0777181053 6362976.949530195, 673355.1857499854 6362983.546437496, 673351.0723593595 6363000, 673593.7562033214 6363000, 673600.5005554566 6362991.191050272))" diff --git a/test/data/donor_infos/donor_info_674_6363_no_donor.csv b/test/data/donor_infos/donor_info_674_6363_no_donor.csv new file mode 100644 index 0000000..f3f4494 --- /dev/null +++ b/test/data/donor_infos/donor_info_674_6363_no_donor.csv @@ -0,0 +1 @@ +,x,y,full_path,geometry diff --git a/test/data/fake_laz_files/data/empty_laz_0673_6363_to_test_pattern_matching_1.laz b/test/data/fake_laz_files/data/empty_laz_0673_6363_to_test_pattern_matching_1.laz new file mode 100644 index 0000000..e69de29 diff --git a/test/data/fake_laz_files/data/empty_laz_0673_6363_to_test_pattern_matching_2.laz b/test/data/fake_laz_files/data/empty_laz_0673_6363_to_test_pattern_matching_2.laz new file mode 100644 index 0000000..e69de29 diff --git a/test/data/lidar_HD_decimated/Semis_2022_0673_6362_LA93_IGN69_decimated.laz b/test/data/lidar_HD_decimated/Semis_2022_0673_6362_LA93_IGN69_decimated.laz new file mode 100644 index 0000000..7edff32 Binary files /dev/null and b/test/data/lidar_HD_decimated/Semis_2022_0673_6362_LA93_IGN69_decimated.laz differ diff --git a/test/data/lidar_HD_decimated/Semis_2022_0673_6363_LA93_IGN69_decimated.laz b/test/data/lidar_HD_decimated/Semis_2022_0673_6363_LA93_IGN69_decimated.laz new file mode 100644 index 0000000..87d54d9 Binary files /dev/null and b/test/data/lidar_HD_decimated/Semis_2022_0673_6363_LA93_IGN69_decimated.laz differ diff --git a/test/data/lidar_HD_decimated/Semis_2022_0674_6363_LA93_IGN69_decimated.laz b/test/data/lidar_HD_decimated/Semis_2022_0674_6363_LA93_IGN69_decimated.laz new file mode 100644 index 0000000..c1147c3 Binary files /dev/null and b/test/data/lidar_HD_decimated/Semis_2022_0674_6363_LA93_IGN69_decimated.laz differ diff --git a/test/data/shapefile_local/patchwork_geometries.dbf b/test/data/shapefile_local/patchwork_geometries.dbf new file mode 100644 index 0000000..e9b8144 Binary files /dev/null and b/test/data/shapefile_local/patchwork_geometries.dbf differ diff --git a/test/data/shapefile_local/patchwork_geometries.shp b/test/data/shapefile_local/patchwork_geometries.shp new file mode 100644 index 0000000..b2f95df Binary files /dev/null and b/test/data/shapefile_local/patchwork_geometries.shp differ diff --git a/test/data/shapefile_local/patchwork_geometries.shx b/test/data/shapefile_local/patchwork_geometries.shx new file mode 100644 index 0000000..59d87a4 Binary files /dev/null and b/test/data/shapefile_local/patchwork_geometries.shx differ diff --git a/test/data/shapefiles_nok/patchwork_geometries_coord_not_in_all_names.dbf b/test/data/shapefiles_nok/patchwork_geometries_coord_not_in_all_names.dbf new file mode 100644 index 0000000..d749a3a Binary files /dev/null and b/test/data/shapefiles_nok/patchwork_geometries_coord_not_in_all_names.dbf differ diff --git a/test/data/shapefiles_nok/patchwork_geometries_coord_not_in_all_names.shp b/test/data/shapefiles_nok/patchwork_geometries_coord_not_in_all_names.shp new file mode 100644 index 0000000..b2f95df Binary files /dev/null and b/test/data/shapefiles_nok/patchwork_geometries_coord_not_in_all_names.shp differ diff --git a/test/data/shapefiles_nok/patchwork_geometries_coord_not_in_all_names.shx b/test/data/shapefiles_nok/patchwork_geometries_coord_not_in_all_names.shx new file mode 100644 index 0000000..59d87a4 Binary files /dev/null and b/test/data/shapefiles_nok/patchwork_geometries_coord_not_in_all_names.shx differ diff --git a/test/data/shapefiles_nok/patchwork_geometries_no_match_for_las_tile.dbf b/test/data/shapefiles_nok/patchwork_geometries_no_match_for_las_tile.dbf new file mode 100644 index 0000000..07701b2 Binary files /dev/null and b/test/data/shapefiles_nok/patchwork_geometries_no_match_for_las_tile.dbf differ diff --git a/test/data/shapefiles_nok/patchwork_geometries_no_match_for_las_tile.shp b/test/data/shapefiles_nok/patchwork_geometries_no_match_for_las_tile.shp new file mode 100644 index 0000000..f15bdf6 Binary files /dev/null and b/test/data/shapefiles_nok/patchwork_geometries_no_match_for_las_tile.shp differ diff --git a/test/data/shapefiles_nok/patchwork_geometries_no_match_for_las_tile.shx b/test/data/shapefiles_nok/patchwork_geometries_no_match_for_las_tile.shx new file mode 100644 index 0000000..f8a458a Binary files /dev/null and b/test/data/shapefiles_nok/patchwork_geometries_no_match_for_las_tile.shx differ diff --git a/test/data/shapefiles_nok/patchwork_geometries_path_does_not_exist.dbf b/test/data/shapefiles_nok/patchwork_geometries_path_does_not_exist.dbf new file mode 100644 index 0000000..5997bd6 Binary files /dev/null and b/test/data/shapefiles_nok/patchwork_geometries_path_does_not_exist.dbf differ diff --git a/test/data/shapefiles_nok/patchwork_geometries_path_does_not_exist.shp b/test/data/shapefiles_nok/patchwork_geometries_path_does_not_exist.shp new file mode 100644 index 0000000..b2f95df Binary files /dev/null and b/test/data/shapefiles_nok/patchwork_geometries_path_does_not_exist.shp differ diff --git a/test/data/shapefiles_nok/patchwork_geometries_path_does_not_exist.shx b/test/data/shapefiles_nok/patchwork_geometries_path_does_not_exist.shx new file mode 100644 index 0000000..59d87a4 Binary files /dev/null and b/test/data/shapefiles_nok/patchwork_geometries_path_does_not_exist.shx differ diff --git a/test/data/shapefiles_nok/patchwork_geometries_several_matches_for_las_tile.dbf b/test/data/shapefiles_nok/patchwork_geometries_several_matches_for_las_tile.dbf new file mode 100644 index 0000000..ef59f6e Binary files /dev/null and b/test/data/shapefiles_nok/patchwork_geometries_several_matches_for_las_tile.dbf differ diff --git a/test/data/shapefiles_nok/patchwork_geometries_several_matches_for_las_tile.shp b/test/data/shapefiles_nok/patchwork_geometries_several_matches_for_las_tile.shp new file mode 100644 index 0000000..f075157 Binary files /dev/null and b/test/data/shapefiles_nok/patchwork_geometries_several_matches_for_las_tile.shp differ diff --git a/test/data/shapefiles_nok/patchwork_geometries_several_matches_for_las_tile.shx b/test/data/shapefiles_nok/patchwork_geometries_several_matches_for_las_tile.shx new file mode 100644 index 0000000..2ca3243 Binary files /dev/null and b/test/data/shapefiles_nok/patchwork_geometries_several_matches_for_las_tile.shx differ diff --git a/test/test_lidar_selecter.py b/test/test_lidar_selecter.py deleted file mode 100644 index dbcb0e8..0000000 --- a/test/test_lidar_selecter.py +++ /dev/null @@ -1,118 +0,0 @@ -import geopandas as gpd -import laspy -import numpy as np -from hydra import compose, initialize -from pandas import DataFrame -from shapely.geometry import MultiPolygon - -import patchwork.constants as c -from patchwork.lidar_selecter import cut_lidar, select_lidar, update_df_result - -CRS = 2154 -TILE_SIZE = 1000 -SHAPE_CORNER_1 = (0, 0) -SHAPE_CORNER_2 = (1000, 0) -SHAPE_CORNER_3 = (1000, 1000) - -POINT_INSIDE_1 = (750, 500, 1) -POINT_INSIDE_2 = (500, 250, 2) -POINT_OUTSIDE_1 = (500, 750, 3) -POINT_OUTSIDE_2 = (250, 500, 4) - -INPUT_DIRECTORY = "las" -OUTPUT_DIRECTORY = "output_directory" -SUBDIRECTORY_NAME = "test" -LASFILE_NAME = "las.las" - - -def test_cut_lidar(): - shapefile_geometry = MultiPolygon( - [ - ([SHAPE_CORNER_1, SHAPE_CORNER_2, SHAPE_CORNER_3],), - ] - ) - las_points = np.array( - [POINT_INSIDE_1, POINT_INSIDE_2, POINT_OUTSIDE_1, POINT_OUTSIDE_2], - dtype=[("x", "float32"), ("y", "float32"), ("z", "float32")], - ) - points_in_geometry = cut_lidar(las_points, shapefile_geometry) - list_points_in_geometry = points_in_geometry.tolist() - assert POINT_INSIDE_1 in list_points_in_geometry - assert POINT_INSIDE_2 in list_points_in_geometry - assert POINT_OUTSIDE_1 not in list_points_in_geometry - assert POINT_OUTSIDE_2 not in list_points_in_geometry - - -def test_select_lidar(tmp_path_factory): - # shapefile creation - shp_dir = tmp_path_factory.mktemp("shapefile") - shp_name = "shapefile.shp" - shapefile_path = shp_dir / shp_name - - shapefile_geometry = MultiPolygon( - [ - ([SHAPE_CORNER_1, SHAPE_CORNER_2, SHAPE_CORNER_3],), - ] - ) - gpd_shapefile_geometry = gpd.GeoDataFrame({"geometry": [shapefile_geometry]}, crs=CRS) - gpd_shapefile_geometry.to_file(shapefile_path) - - # las creation - input_directory = tmp_path_factory.mktemp(INPUT_DIRECTORY) - las_path = input_directory / LASFILE_NAME - - las = laspy.LasData(laspy.LasHeader(point_format=3, version="1.4")) - las_points = np.array( - [POINT_INSIDE_1, POINT_INSIDE_2, POINT_OUTSIDE_1, POINT_OUTSIDE_2], - dtype=[("x", "float32"), ("y", "float32"), ("z", "float32")], - ) - - las.x = las_points["x"] - las.y = las_points["y"] - las.z = las_points["z"] - - las.write(las_path) - - # create output directory - output_directory = las_path = tmp_path_factory.mktemp(OUTPUT_DIRECTORY) - - # create teh dataframe (to put the result in) - data = {c.COORDINATES_KEY: [], c.DONOR_FILE_KEY: [], c.RECIPIENT_FILE_KEY: []} - df_result = DataFrame(data=data) - - with initialize(version_base="1.2", config_path="../configs"): - config = compose( - config_name="configs_patchwork.yaml", - overrides=[f"filepath.SHP_DIRECTORY={shp_dir}", f"filepath.SHP_NAME={shp_name}", f"TILE_SIZE={TILE_SIZE}"], - ) - subdirectory_name = SUBDIRECTORY_NAME - select_lidar(config, input_directory, output_directory, subdirectory_name, df_result, c.DONOR_FILE_KEY, True) - - output_las_path = output_directory / subdirectory_name / LASFILE_NAME - with laspy.open(output_las_path) as las_file: - raw_las_points = las_file.read().points - x = raw_las_points.x - y = raw_las_points.y - z = raw_las_points.z - - las_points_list = list(zip(x, y, z)) - - assert POINT_INSIDE_1 in las_points_list - assert POINT_INSIDE_2 in las_points_list - assert POINT_OUTSIDE_1 not in las_points_list - assert POINT_OUTSIDE_2 not in las_points_list - - -def test_update_df_result(): - data = {c.COORDINATES_KEY: [], c.DONOR_FILE_KEY: [], c.RECIPIENT_FILE_KEY: []} - df_result = DataFrame(data=data) - corner_string = "1111_2222" - donor_path = "dummy_path_1" - df_result = update_df_result(df_result, c.DONOR_FILE_KEY, corner_string, donor_path) - recipient_path = "dummy_path_2" - df_result = update_df_result(df_result, c.RECIPIENT_FILE_KEY, corner_string, recipient_path) - assert len(df_result) == 1 - row_0 = df_result.loc[0] - assert row_0[c.COORDINATES_KEY] == corner_string - assert row_0[c.DONOR_FILE_KEY] == donor_path - assert row_0[c.RECIPIENT_FILE_KEY] == recipient_path diff --git a/test/test_patchwork.py b/test/test_patchwork.py index 30aad4f..5623af6 100644 --- a/test/test_patchwork.py +++ b/test/test_patchwork.py @@ -1,24 +1,22 @@ import os +import geopandas as gpd import laspy import numpy as np import pandas as pd import pytest from hydra import compose, initialize -from pandas import DataFrame +from pdaltools.las_info import get_tile_origin_using_header_info import patchwork.constants as c from patchwork.patchwork import ( append_points, get_complementary_points, - get_donor_from_csv, - get_donor_path, get_field_from_header, get_selected_classes_points, get_type, patchwork, ) -from patchwork.tools import get_tile_origin_from_pointcloud RECIPIENT_TEST_DIR = "test/data/" RECIPIENT_TEST_NAME = "recipient_test.laz" @@ -32,20 +30,7 @@ NEW_COLUMN_SIZE = 8 VALUE_ADDED_POINTS = 1 -DONOR_TEST_DIR = "test/data/" -DONOR_TEST_NAME = "donor_test.las" - -RECIPIENT_MORE_FIELDS_TEST_DIR = "test/data" -RECIPIENT_MORE_FIELDS_TEST_NAME = "recipient_more_fields_test.laz" - -DONOR_MORE_FIELDS_TEST_DIR = "test/data" -DONOR_MORE_FIELDS_TEST_NAME = "donor_more_fields_test.las" - -RECIPIENT_SLIDED_TEST_DIR = "test/data" -RECIPIENT_SLIDED_TEST_NAME = "recipient_slided_test.laz" - - -COORDINATES = "1234_6789" +SHP_X_Y_TO_METER_FACTOR = 1000 def test_get_field_from_header(): @@ -57,7 +42,6 @@ def test_get_field_from_header(): def test_get_selected_classes_points(): - with initialize(version_base="1.2", config_path="../configs"): config = compose( config_name="configs_patchwork.yaml", @@ -67,14 +51,10 @@ def test_get_selected_classes_points(): f"RECIPIENT_CLASS_LIST={RECIPIENT_CLASS_LIST}", ], ) - - with laspy.open( - os.path.join(config.filepath.RECIPIENT_DIRECTORY, config.filepath.RECIPIENT_NAME) - ) as recipient_file: + recipient_path = os.path.join(config.filepath.RECIPIENT_DIRECTORY, config.filepath.RECIPIENT_NAME) + tile_origin_recipient = get_tile_origin_using_header_info(recipient_path, config.TILE_SIZE) + with laspy.open(recipient_path) as recipient_file: recipient_points = recipient_file.read().points - - tile_origin_recipient = get_tile_origin_from_pointcloud(config, recipient_points) - df_recipient_points = get_selected_classes_points( config, tile_origin_recipient, recipient_points, config.RECIPIENT_CLASS_LIST, [] ) @@ -82,82 +62,113 @@ def test_get_selected_classes_points(): assert classification in RECIPIENT_CLASS_LIST -def test_get_complementary_points(): +@pytest.mark.parametrize( + "donor_info_path, recipient_path, x, y, expected_nb_points", + # expected_nb_points value set after inspection of the initial result using qgis: + # - there are points only inside the shapefile geometry + # - when visualizing a grid, there seems to be no points in the cells where there is ground points in the + # recipient laz + [ + ( + "test/data/donor_infos/donor_info_673_6362_one_donor.csv", + "test/data/lidar_HD_decimated/Semis_2022_0673_6362_LA93_IGN69_decimated.laz", + 673, + 6362, + 128675, + ), + ( + "test/data/donor_infos/donor_info_673_6363_two_donors.csv", + "test/data/lidar_HD_decimated/Semis_2022_0673_6363_LA93_IGN69_decimated.laz", + 673, + 6363, + 149490, + ), + ( + "test/data/donor_infos/donor_info_674_6363_no_donor.csv", + "test/data/lidar_HD_decimated/Semis_2022_0674_6363_LA93_IGN69_decimated.laz", + 674, + 6363, + 0, + ), + ], +) +def test_get_complementary_points(donor_info_path, recipient_path, x, y, expected_nb_points): + df = pd.read_csv(donor_info_path, encoding="utf-8") + s = gpd.GeoSeries.from_wkt(df.geometry) + df_donor_info = gpd.GeoDataFrame(data=df, geometry=s) + with initialize(version_base="1.2", config_path="../configs"): config = compose( config_name="configs_patchwork.yaml", overrides=[ - f"filepath.DONOR_DIRECTORY={DONOR_TEST_DIR}", - f"filepath.DONOR_NAME={DONOR_TEST_NAME}", - f"filepath.RECIPIENT_DIRECTORY={RECIPIENT_TEST_DIR}", - f"filepath.RECIPIENT_NAME={RECIPIENT_TEST_NAME}", f"DONOR_CLASS_LIST={DONOR_CLASS_LIST}", f"RECIPIENT_CLASS_LIST={RECIPIENT_CLASS_LIST}", f"+VIRTUAL_CLASS_TRANSLATION={VIRTUAL_CLASS_TRANSLATION}", ], ) + complementary_points = get_complementary_points(df_donor_info, recipient_path, (x, y), config) + + assert np.all(complementary_points["x"] >= x * SHP_X_Y_TO_METER_FACTOR) + assert np.all(complementary_points["x"] <= (x + 1) * SHP_X_Y_TO_METER_FACTOR) + assert np.all(complementary_points["y"] >= (y - 1) * SHP_X_Y_TO_METER_FACTOR) + assert np.all(complementary_points["y"] <= y * SHP_X_Y_TO_METER_FACTOR) - complementary_points = get_complementary_points(config) - assert len(complementary_points) == 320 + assert len(complementary_points.index) == expected_nb_points -def test_get_complementary_points_2(): +def test_get_complementary_points_2_more_fields(tmp_path_factory): """test selected_classes_points with more fields in files, different from each other's""" + original_recipient_path = "test/data/lidar_HD_decimated/Semis_2022_0673_6362_LA93_IGN69_decimated.laz" + original_donor_path = ( + "test/data/aveyron_aval_lidarBD/data/NUALID_1-0_VLIDAVEYRONAVAL_PTS_0673_6362_LAMB93_IGN69_20210319.laz" + ) + original_donor_info_path = "test/data/donor_infos/donor_info_673_6362_one_donor.csv" + x = 673 + y = 6362 + + tmp_dir = tmp_path_factory.mktemp("data") + tmp_recipient_name = "recipient_with_extra_dims.laz" + tmp_recipient_path = os.path.join(tmp_dir, tmp_recipient_name) + + tmp_donor_name = "donor_with_extra_dims.laz" + tmp_donor_path = os.path.join(tmp_dir, tmp_donor_name) + + df = pd.read_csv(original_donor_info_path, encoding="utf-8") + s = gpd.GeoSeries.from_wkt(df.geometry) + df_donor_info = gpd.GeoDataFrame(data=df, geometry=s) + df_donor_info["full_path"] = [tmp_donor_path] + extra_fields_for_recipient = ["f1", "f2"] - las = laspy.read(os.path.join(RECIPIENT_TEST_DIR, RECIPIENT_TEST_NAME)) + las = laspy.read(original_recipient_path) for field in extra_fields_for_recipient: las.add_extra_dim(laspy.ExtraBytesParams(name=field, type=np.uint64)) - las.write(os.path.join(RECIPIENT_MORE_FIELDS_TEST_DIR, RECIPIENT_MORE_FIELDS_TEST_NAME)) + las.write(tmp_recipient_path) extra_fields_for_donor = ["f3", "f4"] - las = laspy.read(os.path.join(DONOR_TEST_DIR, DONOR_TEST_NAME)) + las = laspy.read(original_donor_path) for field in extra_fields_for_donor: las.add_extra_dim(laspy.ExtraBytesParams(name=field, type=np.uint64)) - las.write(os.path.join(DONOR_MORE_FIELDS_TEST_DIR, DONOR_MORE_FIELDS_TEST_NAME)) + las.write(tmp_donor_path) with initialize(version_base="1.2", config_path="../configs"): config = compose( config_name="configs_patchwork.yaml", overrides=[ - f"filepath.RECIPIENT_DIRECTORY={RECIPIENT_MORE_FIELDS_TEST_DIR}", - f"filepath.RECIPIENT_NAME={RECIPIENT_MORE_FIELDS_TEST_NAME}", - f"filepath.DONOR_DIRECTORY={DONOR_MORE_FIELDS_TEST_DIR}", - f"filepath.DONOR_NAME={DONOR_MORE_FIELDS_TEST_NAME}", f"DONOR_CLASS_LIST={DONOR_CLASS_LIST}", f"RECIPIENT_CLASS_LIST={RECIPIENT_CLASS_LIST}", f"+VIRTUAL_CLASS_TRANSLATION={VIRTUAL_CLASS_TRANSLATION}", ], ) - complementary_points = get_complementary_points(config) - assert len(complementary_points) == 320 - columns = complementary_points.columns - for field in extra_fields_for_recipient: # no extra field from the recipient should exist - assert field not in columns - for field in extra_fields_for_donor: # every extra field from the donor should exist... - assert field in columns - assert complementary_points[field].all() == 0 # ...but should be at 0 - - -def test_get_complementary_points_3(): - """test selected_classes_points with 2 files from different areas""" - with initialize(version_base="1.2", config_path="../configs"): - config = compose( - config_name="configs_patchwork.yaml", - overrides=[ - f"filepath.DONOR_DIRECTORY={DONOR_TEST_DIR}", - f"filepath.DONOR_NAME={DONOR_TEST_NAME}", - f"filepath.RECIPIENT_DIRECTORY={RECIPIENT_SLIDED_TEST_DIR}", - f"filepath.RECIPIENT_NAME={RECIPIENT_SLIDED_TEST_NAME}", - ], - ) - - las = laspy.read(os.path.join(RECIPIENT_TEST_DIR, RECIPIENT_TEST_NAME)) - las.points["x"] = las.points["x"] + config.TILE_SIZE - las.write(os.path.join(RECIPIENT_SLIDED_TEST_DIR, RECIPIENT_SLIDED_TEST_NAME)) + complementary_points = get_complementary_points(df_donor_info, tmp_recipient_path, (x, y), config) - with pytest.raises(Exception): - get_complementary_points(config) + assert len(complementary_points.index) == 128675 + columns = complementary_points.columns + for field in extra_fields_for_recipient: # no extra field from the recipient should exist + assert field not in columns + for field in extra_fields_for_donor: # every extra field from the donor should exist... + assert field in columns + assert complementary_points[field].all() == 0 # ...but should be at 0 def test_get_type(): @@ -273,160 +284,125 @@ def test_append_points_new_column(tmp_path_factory): assert max(new_column[:-2]) == 0 -def test_get_donor_from_csv(tmp_path_factory): - csv_file_path = tmp_path_factory.mktemp("csv") / "recipients_donors_links.csv" - donor_more_fields_test_path = os.path.join(DONOR_MORE_FIELDS_TEST_DIR, DONOR_MORE_FIELDS_TEST_NAME) - recipient_more_fields_test_path = os.path.join(RECIPIENT_TEST_DIR, RECIPIENT_TEST_NAME) - data = { - c.COORDINATES_KEY: [ - COORDINATES, - ], - c.DONOR_FILE_KEY: [ - donor_more_fields_test_path, - ], - c.RECIPIENT_FILE_KEY: [ - recipient_more_fields_test_path, - ], - } - DataFrame(data=data).to_csv(csv_file_path) - - donor_file_path = get_donor_from_csv(recipient_more_fields_test_path, csv_file_path) - assert donor_file_path == donor_more_fields_test_path - - -def test_get_donor_path(tmp_path_factory): - # check get_donor_path when no csv - with initialize(version_base="1.2", config_path="../configs"): - config = compose( - config_name="configs_patchwork.yaml", - overrides=[ - f"filepath.DONOR_DIRECTORY={DONOR_TEST_DIR}", - f"filepath.DONOR_NAME={DONOR_TEST_NAME}", - f"filepath.RECIPIENT_DIRECTORY={RECIPIENT_SLIDED_TEST_DIR}", - f"filepath.RECIPIENT_NAME={RECIPIENT_SLIDED_TEST_NAME}", - ], - ) - donor_dir, donor_name = get_donor_path(config) - assert donor_dir == DONOR_TEST_DIR - assert donor_name == DONOR_TEST_NAME - - # check get_donor_path when csv but no matching donor in it - csv_file_dir = tmp_path_factory.mktemp("csv") - csv_file_name = "recipients_donors_links.csv" - csv_file_path = os.path.join(csv_file_dir, csv_file_name) - - data = {c.COORDINATES_KEY: [], c.DONOR_FILE_KEY: [], c.RECIPIENT_FILE_KEY: []} - DataFrame(data=data).to_csv(csv_file_path) - - with initialize(version_base="1.2", config_path="../configs"): - config = compose( - config_name="configs_patchwork.yaml", - overrides=[ - f"filepath.DONOR_DIRECTORY={DONOR_TEST_DIR}", - f"filepath.DONOR_NAME={DONOR_TEST_NAME}", - f"filepath.RECIPIENT_DIRECTORY={RECIPIENT_SLIDED_TEST_DIR}", - f"filepath.RECIPIENT_NAME={RECIPIENT_SLIDED_TEST_NAME}", - f"filepath.CSV_DIRECTORY={csv_file_dir}", - f"filepath.CSV_NAME={csv_file_name}", - ], - ) - - donor_dir, donor_name = get_donor_path(config) - assert donor_dir == "" - assert donor_name == "" - - # check get_donor_path when csv but with a matching donor in it - donor_more_fields_test_path = os.path.join(DONOR_MORE_FIELDS_TEST_DIR, DONOR_MORE_FIELDS_TEST_NAME) - recipient_more_fields_test_path = os.path.join(RECIPIENT_TEST_DIR, RECIPIENT_TEST_NAME) - data = { - c.COORDINATES_KEY: [ - COORDINATES, - ], - c.DONOR_FILE_KEY: [ - donor_more_fields_test_path, - ], - c.RECIPIENT_FILE_KEY: [ - recipient_more_fields_test_path, - ], - } - DataFrame(data=data).to_csv(csv_file_path) - - with initialize(version_base="1.2", config_path="../configs"): - config = compose( - config_name="configs_patchwork.yaml", - overrides=[ - f"filepath.DONOR_DIRECTORY={DONOR_TEST_DIR}", - f"filepath.DONOR_NAME={DONOR_TEST_NAME}", - f"filepath.RECIPIENT_DIRECTORY={RECIPIENT_TEST_DIR}", - f"filepath.RECIPIENT_NAME={RECIPIENT_TEST_NAME}", - f"filepath.CSV_DIRECTORY={csv_file_dir}", - f"filepath.CSV_NAME={csv_file_name}", - ], - ) - - donor_dir, donor_name = get_donor_path(config) - assert donor_dir == DONOR_MORE_FIELDS_TEST_DIR - assert donor_name == DONOR_MORE_FIELDS_TEST_NAME - - -def test_patchwork_default(tmp_path_factory): +@pytest.mark.parametrize( + "recipient_path, expected_nb_added_points", + # expected_nb_points value set after inspection of the initial result using qgis: + # - there are points only inside the shapefile geometry + # - when visualizing a grid, there seems to be no points in the cells where there is ground points in the + # recipient laz + [ + ( + "test/data/lidar_HD_decimated/Semis_2022_0673_6362_LA93_IGN69_decimated.laz", + 128675, + ), # One donor + ( + "test/data/lidar_HD_decimated/Semis_2022_0673_6363_LA93_IGN69_decimated.laz", + 149490, + ), # Two donors + ( + "test/data/lidar_HD_decimated/Semis_2022_0674_6363_LA93_IGN69_decimated.laz", + 0, + ), # No donor + ], +) +def test_patchwork_default(tmp_path_factory, recipient_path, expected_nb_added_points): + input_shp_path = "test/data/shapefile_local/patchwork_geometries.shp" tmp_file_dir = tmp_path_factory.mktemp("data") tmp_output_las_name = "result_patchwork.laz" - tmp_output_indices_map_name = "result_patchwerk_indices.tif" + tmp_output_indices_map_name = "result_patchwork_indices.tif" with initialize(version_base="1.2", config_path="../configs"): config = compose( config_name="configs_patchwork.yaml", overrides=[ - f"filepath.RECIPIENT_DIRECTORY={RECIPIENT_TEST_DIR}", - f"filepath.RECIPIENT_NAME={RECIPIENT_TEST_NAME}", - f"filepath.DONOR_DIRECTORY={DONOR_TEST_DIR}", - f"filepath.DONOR_NAME={DONOR_TEST_NAME}", + f"filepath.RECIPIENT_DIRECTORY={os.path.dirname(recipient_path)}", + f"filepath.RECIPIENT_NAME={os.path.basename(recipient_path)}", + f"filepath.SHP_DIRECTORY={os.path.dirname(input_shp_path)}", + f"filepath.SHP_NAME={os.path.basename(input_shp_path)}", f"filepath.OUTPUT_DIR={tmp_file_dir}", f"filepath.OUTPUT_NAME={tmp_output_las_name}", f"filepath.OUTPUT_INDICES_MAP_DIR={tmp_file_dir}", f"filepath.OUTPUT_INDICES_MAP_NAME={tmp_output_indices_map_name}", + f"DONOR_CLASS_LIST={DONOR_CLASS_LIST}", + f"RECIPIENT_CLASS_LIST={RECIPIENT_CLASS_LIST}", + f"+VIRTUAL_CLASS_TRANSLATION={VIRTUAL_CLASS_TRANSLATION}", + "NEW_COLUMN=null", ], ) patchwork(config) - recipient_path = os.path.join(config.filepath.RECIPIENT_DIRECTORY, config.filepath.RECIPIENT_NAME) output_path = os.path.join(tmp_file_dir, tmp_output_las_name) indices_map_path = os.path.join(tmp_file_dir, tmp_output_indices_map_name) assert os.path.isfile(output_path) assert os.path.isfile(indices_map_path) + with laspy.open(recipient_path) as las_file: recipient_points = las_file.read().points with laspy.open(output_path) as las_file: output_points = las_file.read().points - assert len(output_points) > len(recipient_points) - - -def test_patchwork_empty_donor(tmp_path_factory): + assert {n for n in las_file.header.point_format.dimension_names} == { + n for n in las_file.header.point_format.standard_dimension_names + } + assert len(output_points) == len(recipient_points) + expected_nb_added_points + + +@pytest.mark.parametrize( + "recipient_path, expected_nb_added_points", + # expected_nb_points value set after inspection of the initial result using qgis: + # - there are points only inside the shapefile geometry + # - when visualizing a grid, there seems to be no points in the cells where there is ground points in the + # recipient laz + [ + ( + "test/data/lidar_HD_decimated/Semis_2022_0673_6362_LA93_IGN69_decimated.laz", + 128675, + ), # One donor + ( + "test/data/lidar_HD_decimated/Semis_2022_0673_6363_LA93_IGN69_decimated.laz", + 149490, + ), # Two donors + ( + "test/data/lidar_HD_decimated/Semis_2022_0674_6363_LA93_IGN69_decimated.laz", + 0, + ), # No donor + ], +) +def test_patchwork_with_origin(tmp_path_factory, recipient_path, expected_nb_added_points): + input_shp_path = "test/data/shapefile_local/patchwork_geometries.shp" tmp_file_dir = tmp_path_factory.mktemp("data") tmp_output_las_name = "result_patchwork.laz" - tmp_output_indices_map_name = "result_patchwerk_indices.tif" + tmp_output_indices_map_name = "result_patchwork_indices.tif" with initialize(version_base="1.2", config_path="../configs"): config = compose( config_name="configs_patchwork.yaml", overrides=[ - f"filepath.RECIPIENT_DIRECTORY={RECIPIENT_TEST_DIR}", - f"filepath.RECIPIENT_NAME={RECIPIENT_TEST_NAME}", + f"filepath.RECIPIENT_DIRECTORY={os.path.dirname(recipient_path)}", + f"filepath.RECIPIENT_NAME={os.path.basename(recipient_path)}", + f"filepath.SHP_DIRECTORY={os.path.dirname(input_shp_path)}", + f"filepath.SHP_NAME={os.path.basename(input_shp_path)}", f"filepath.OUTPUT_DIR={tmp_file_dir}", f"filepath.OUTPUT_NAME={tmp_output_las_name}", f"filepath.OUTPUT_INDICES_MAP_DIR={tmp_file_dir}", f"filepath.OUTPUT_INDICES_MAP_NAME={tmp_output_indices_map_name}", + f"DONOR_CLASS_LIST={DONOR_CLASS_LIST}", + f"RECIPIENT_CLASS_LIST={RECIPIENT_CLASS_LIST}", + "NEW_COLUMN='Origin'", ], ) patchwork(config) - recipient_path = os.path.join(config.filepath.RECIPIENT_DIRECTORY, config.filepath.RECIPIENT_NAME) output_path = os.path.join(tmp_file_dir, tmp_output_las_name) indices_map_path = os.path.join(tmp_file_dir, tmp_output_indices_map_name) assert os.path.isfile(output_path) assert os.path.isfile(indices_map_path) + with laspy.open(recipient_path) as las_file: recipient_points = las_file.read().points with laspy.open(output_path) as las_file: output_points = las_file.read().points + assert {n for n in las_file.header.point_format.dimension_names} == { + n for n in las_file.header.point_format.standard_dimension_names + } | {"Origin"} - assert len(output_points) == len(recipient_points) + assert len(output_points) == len(recipient_points) + expected_nb_added_points + assert np.sum(output_points.Origin == 0) == len(recipient_points) + assert np.sum(output_points.Origin == 1) == expected_nb_added_points diff --git a/test/test_shapefile_data_extraction.py b/test/test_shapefile_data_extraction.py new file mode 100644 index 0000000..977e4c6 --- /dev/null +++ b/test/test_shapefile_data_extraction.py @@ -0,0 +1,67 @@ +import pytest + +from patchwork.shapefile_data_extraction import get_donor_info_from_shapefile + +INPUT_SHP_PATH = "test/data/shapefile_local/patchwork_geometries.shp" +DONOR_SUBDIRECTORY = "data" + + +@pytest.mark.parametrize( + "x,y, expected_full_path", + [ + ( + 673, + 6362, + {"test/data/aveyron_aval_lidarBD/data/NUALID_1-0_VLIDAVEYRONAVAL_PTS_0673_6362_LAMB93_IGN69_20210319.laz"}, + ), # Expect only one output file + ( + 673, + 6363, + { + "test/data/aveyron_aval_lidarBD/data/NUALID_1-0_VLIDAVEYRONAVAL_PTS_0673_6363_LAMB93_IGN69_20210319.laz", # noqa: E501 + "test/data/aveyron_lidarBD/data/NUALID_1-0_IAVEY_PTS_0673_6363_LAMB93_IGN69_20170519.laz", + }, + ), # Expect 2 output files + (673, 6365, set()), # Expect no output file + ], +) +def test_get_donor_info_from_shapefile(x, y, expected_full_path): + gdf = get_donor_info_from_shapefile(INPUT_SHP_PATH, x, y, DONOR_SUBDIRECTORY) + # Check that all paths are filled + assert set(gdf.columns) == {"x", "y", "full_path", "geometry"} + assert len(gdf.index) == len(expected_full_path) + assert set(gdf["full_path"]) == expected_full_path + + +@pytest.mark.parametrize( + "input_shp, x, y, error_type", + [ + ( + "test/data/shapefiles_nok/patchwork_geometries_coord_not_in_all_names.shp", + 673, + 6362, + NotImplementedError, + ), + ( + "test/data/shapefiles_nok/patchwork_geometries_path_does_not_exist.shp", + 673, + 6363, + FileNotFoundError, + ), + ( + "test/data/shapefiles_nok/patchwork_geometries_no_match_for_las_tile.shp", + 673, + 6365, + FileNotFoundError, + ), # Tile described in shapefile but not found in directory + ( + "test/data/shapefiles_nok/patchwork_geometries_several_matches_for_las_tile.shp", + 673, + 6363, + RuntimeError, + ), + ], +) +def test_get_donor_info_from_shapefile_raise_error(input_shp, x, y, error_type): + with pytest.raises(error_type): + get_donor_info_from_shapefile(input_shp, x, y, DONOR_SUBDIRECTORY) diff --git a/test/test_tools.py b/test/test_tools.py deleted file mode 100644 index 58514f6..0000000 --- a/test/test_tools.py +++ /dev/null @@ -1,61 +0,0 @@ -import numpy as np -import pytest -from hydra import compose, initialize - -from patchwork.tools import get_tile_origin_from_pointcloud - -TILE_SIZE = 1000 - - -def test_get_tile_origin_from_pointcloud(): - - with initialize(version_base="1.2", config_path="../configs"): - config = compose( - config_name="configs_patchwork.yaml", - overrides=[ - f"TILE_SIZE={TILE_SIZE}", - ], - ) - - # basic test - list_x = [1100, 1500] - list_y = [2200, 2800] - points = np.core.records.fromarrays([list_x, list_y], names="x,y") - - corner_x, corner_y = get_tile_origin_from_pointcloud(config, points) - assert corner_x == 1000 - assert corner_y == 3000 - - # limit test 1 - list_x = [1000, 2000] - list_y = [1000, 2000] - points = np.core.records.fromarrays([list_x, list_y], names="x,y") - - corner_x, corner_y = get_tile_origin_from_pointcloud(config, points) - assert corner_x == 1000 - assert corner_y == 2000 - - # limit test 2 - list_x = [1500] - list_y = [2300] - points = np.core.records.fromarrays([list_x, list_y], names="x,y") - - corner_x, corner_y = get_tile_origin_from_pointcloud(config, points) - assert corner_x == 1000 - assert corner_y == 3000 - - # limit test 3 - list_x = [] - list_y = [] - points = np.core.records.fromarrays([list_x, list_y], names="x,y") - - with pytest.raises(ValueError): - get_tile_origin_from_pointcloud(config, points) - - # failed test - list_x = [1100, 1500] - list_y = [2200, 3800] - points = np.core.records.fromarrays([list_x, list_y], names="x,y") - - with pytest.raises(ValueError): - get_tile_origin_from_pointcloud(config, points)