diff --git a/.gitignore b/.gitignore index 99b7194..7414c1a 100644 --- a/.gitignore +++ b/.gitignore @@ -4,6 +4,7 @@ __pycache__ test/tmp/ *.log +*.copc.laz # Library deployment *.egg-info diff --git a/CHANGELOG.md b/CHANGELOG.md index 7b61476..5f876ae 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,12 +1,13 @@ # dev -use pdal>=2.9 and deactivate the Dockerfile.pdal from CI +- use pdal>=2.9 and deactivate the Dockerfile.pdal from CI +- [BREAKING CHANGE] Add a [main.py](las_digital_models/main.py) to run the whole pipeline at once (buffer, DTM, DSM, DHM). Use temporary folders for intermediate values. Refactor config files. # v2.1.1 -fix sur le déploiement de l'image Docker +fix CI to deploy docker image # v2.1.0 Custom PDAL: in the docker image, compile custom PDAL (waiting for PDAL 2.9) -fix run_extract_z_virtual_lines_from_raster: output geometries are only LineString (no more MultiLineString) +fix run_extract_z_virtual_lines_from_raster: output geometries are only LineString (no more MultiLineString) # v2.0.0 Rename produit_derive_lidar to las_digital_models diff --git a/README.md b/README.md index 50456e7..57d7770 100644 --- a/README.md +++ b/README.md @@ -124,7 +124,7 @@ the `configs` folder. > * `tile_geometry.tile_width` must contain the tile size in meters > * `tile_geometry.tile_coord_scale` must contain the `coord_x` and `coord_y` scale so that `coord_{x or y} * tile_geometry.tile_coord_scale` are the coordinates of the top-left corner in meters -### Whole pipeline +### Whole pipeline on a folder To run the whole pipeline (DSM + DTM + DHM) on all the LAS files in a folder, use `run.sh`. @@ -147,6 +147,24 @@ It will generate: * ${OUTPUT_DIR}/DSM * ${OUTPUT_DIR}/DHM +### Whole pipeline on a single file + +To run the whole pipeline (buffer + DTM + DSM + DHM) on a single file: + +```bash +python -m las_digital_models.main \ + io.input_dir=INPUT_DIR \ + io.input_filename=INPUT_FILENAME \ + io.output_dir=OUTPUT_DIR \ + tile_geometry.pixel_size=${PIXEL_SIZE} + buffer.size=10 +``` +Any of DTM, DSM or DHM computation can be deactivated using `tasks.dtm=false`, `tasks.dsm=false` +or `tasks.dhm=false`. + +Any other parameter in the `./configs` tree can be overriden in the command (see the doc of +[hydra](https://hydra.cc/) for more details on usage) + ### Buffer To add a buffer to a point cloud using `ign-pdal-tools`: @@ -155,8 +173,7 @@ To add a buffer to a point cloud using `ign-pdal-tools`: python -m las_digital_models.filter_one_tile \ io.input_dir=INPUT_DIR \ io.input_filename=INPUT_FILENAME \ - io.output_dir=OUTPUT_DIR \ - buffer.size=10 + io.output_dir=OUTPUT_DIR ``` Any other parameter in the `./configs` tree can be overriden in the command (see the doc of @@ -172,8 +189,8 @@ python -m las_digital_models.ip_one_tile \ io.input_filename={} \ io.output_dir=${DXM_DIR} \ tile_geometry.pixel_size=${PIXEL_SIZE} \ - filter.dimension="Classification" \ - filter.keep_values=[2,66] + interpolation.custom.filter.dimension="Classification" \ + interpolation.custom.filter.keep_values=[2,66] ``` `filter.keep_values` must be a list inside `[]`, separated by `,` without spaces. diff --git a/configs/buffer/default.yaml b/configs/buffer/default.yaml index 2914e2e..992ffd2 100644 --- a/configs/buffer/default.yaml +++ b/configs/buffer/default.yaml @@ -1 +1,3 @@ -size: 100 # bufer size in meters \ No newline at end of file +size: 100 # buffer size in meters + +output_subdir: null # subdirectory of oi.output_dir in which to save files with buffer, for debug only \ No newline at end of file diff --git a/configs/buffer/test.yaml b/configs/buffer/test.yaml index a316d36..6e0bd6d 100644 --- a/configs/buffer/test.yaml +++ b/configs/buffer/test.yaml @@ -1 +1,3 @@ -size: 10 \ No newline at end of file +size: 10 + +output_subdir: null # subdirectory of oi.output_dir in which to save files with buffer, for debug only \ No newline at end of file diff --git a/configs/config.yaml b/configs/config.yaml index 6b252f6..ccb0eb5 100644 --- a/configs/config.yaml +++ b/configs/config.yaml @@ -6,15 +6,19 @@ # learn more here: https://hydra.cc/docs/next/tutorials/basic/running_your_app/working_directory work_dir: ${hydra:runtime.cwd} +tasks: + dtm: true + dsm: true + dhm: true + # specify here default training configuration defaults: - _self_ # for hydra legacy reasons - - buffer: default.yaml - - filter: dtm.yaml + - interpolation: default.yaml + - dhm: default.yaml - io: default.yaml - tile_geometry: default.yaml # describes input features and classes - - dhm: default.yaml - extract_stat: default.yaml # disable hydra logging diff --git a/configs/dhm/default.yaml b/configs/dhm/default.yaml index a022f28..90d8497 100644 --- a/configs/dhm/default.yaml +++ b/configs/dhm/default.yaml @@ -1,2 +1,9 @@ -input_dsm_dir: /path/to/dsm/dir -input_dtm_dir: /path/to/dtm/dir \ No newline at end of file +output_subfolder: "DHM" + +# To be used only when using dhm_one_tile directly (standalone mode) to run +# DHM from any input DTM and DSM +# DSM and DTM are associated to the tile by their names, which should be like: +# f"{input filename without extension}_{pixel size with unit}" +input_dtm_dir: /path/to/input/dtm/folder +input_dsm_dir: /path/to/input/dsm/folder + diff --git a/configs/dhm/test.yaml b/configs/dhm/test.yaml index d2e15f1..027cbee 100644 --- a/configs/dhm/test.yaml +++ b/configs/dhm/test.yaml @@ -1,2 +1,9 @@ +output_subfolder: "DHM" + +# To be used only when using dhm_one_tile directly (standalone mode) to run +# DHM from any input DTM and DSM +# DSM and DTM are associated to the tile by their names, which should be like: +# f"{input filename without extension}_{pixel size with unit}" +input_dtm_dir: ./test/data/DTM input_dsm_dir: ./test/data/DSM -input_dtm_dir: ./test/data/DTM \ No newline at end of file + diff --git a/configs/filter/default_test.yaml b/configs/filter/default_test.yaml deleted file mode 100644 index 7ac1630..0000000 --- a/configs/filter/default_test.yaml +++ /dev/null @@ -1,2 +0,0 @@ -dimension: Classification -keep_values: [2, 66] # ground + virtual points \ No newline at end of file diff --git a/configs/filter/dsm.yaml b/configs/filter/dsm.yaml deleted file mode 100644 index 5e48120..0000000 --- a/configs/filter/dsm.yaml +++ /dev/null @@ -1,2 +0,0 @@ -dimension: Classification -keep_values: [2, 3, 4, 5, 6, 9, 17] # Classes to keep when running DSM \ No newline at end of file diff --git a/configs/filter/dtm.yaml b/configs/filter/dtm.yaml deleted file mode 100644 index cadb701..0000000 --- a/configs/filter/dtm.yaml +++ /dev/null @@ -1,2 +0,0 @@ -dimension: Classification -keep_values: [2, 9, 66] # Classes to keep when running DTM \ No newline at end of file diff --git a/configs/interpolation/default.yaml b/configs/interpolation/default.yaml new file mode 100644 index 0000000..d334875 --- /dev/null +++ b/configs/interpolation/default.yaml @@ -0,0 +1,19 @@ +dsm: + filter: + dimension: Classification + keep_values: [2, 3, 4, 5, 6, 9, 17] # Classes to keep when running DSM + output_subfolder: "DSM" + +dtm: + filter: + dimension: Classification + keep_values: [2, 9, 66] # Classes to keep when running DTM + output_subfolder: "DTM" + + +# To be used only when using ip_one_tile directly (standalone mode) to run +# interpolation on any kind of filtered data +custom: + filter: + dimension: Classification + keep_values: [2, 9, 66] # Classes to keep when running ip_one_tile diff --git a/configs/interpolation/test.yaml b/configs/interpolation/test.yaml new file mode 100644 index 0000000..78afffd --- /dev/null +++ b/configs/interpolation/test.yaml @@ -0,0 +1,19 @@ +dsm: + filter: + dimension: Classification + keep_values: [2, 3, 4, 5, 6, 9, 17] # Classes to keep when running DSM + output_subfolder: "DSM" + +dtm: + filter: + dimension: Classification + keep_values: [2, 9, 66] # Classes to keep when running DTM + output_subfolder: "DTM" + + +# To be used only when using ip_one_tile directly (standalone mode) to run +# interpolation on any kind of filtered data +custom: + filter: + dimension: Classification + keep_values: [2, 66] # ground + virtual points \ No newline at end of file diff --git a/configs/io/default.yaml b/configs/io/default.yaml index 1d1a6c6..fb19bf3 100644 --- a/configs/io/default.yaml +++ b/configs/io/default.yaml @@ -6,13 +6,14 @@ input_filename: filename.las # Shapefile used to burn no-data value (inside its polygons) no_data_mask_shapefile: null # /path/to/mask.shp -# Force input/output files extension to .las or .laz. If None or not set, use original extension -# Used for output for filter -# Used for input + output in buffer -# Use for input in interpolation -forced_intermediate_ext: null # can be "las", "laz" or null. (intermediate results extension) - # Spatial reference to use to override the one from input las. spatial_reference: EPSG:2154 -output_dir: /path/to/output/folder # Directory folder for saving the outputs \ No newline at end of file +output_dir: /path/to/output/folder # Directory folder for saving the outputs + +# To be used only when using dhm_one_tile directly (standalone mode) to run +# DHM from any input DTM and DSM +# DSM and DTM are associated to the tile by their names, which should be like: +# f"{input filename without extension}_{pixel size with unit}" +input_dtm_dir: /path/to/input/dtm/folder +input_dsm_dir: /path/to/input/dsm/folder \ No newline at end of file diff --git a/configs/io/test.yaml b/configs/io/test.yaml index 0ffda7b..24654c9 100644 --- a/configs/io/test.yaml +++ b/configs/io/test.yaml @@ -6,4 +6,11 @@ forced_intermediate_ext: null spatial_reference: EPSG:2154 -output_dir: ./test/tmp/hydra_ip \ No newline at end of file +output_dir: ./test/tmp/hydra_ip + +# To be used only when using dhm_one_tile directly (standalone mode) to run +# DHM from any input DTM and DSM +# DSM and DTM are associated to the tile by their names, which should be like: +# f"{input filename without extension}_{pixel size with unit}" +input_dsm_dir: ./test/data/DSM +input_dtm_dir: ./test/data/DTM \ No newline at end of file diff --git a/configs/test.yaml b/configs/test.yaml index df87c90..e8f1524 100644 --- a/configs/test.yaml +++ b/configs/test.yaml @@ -3,14 +3,19 @@ # Special config used for tests (test_run_script.py) work_dir: ${hydra:runtime.cwd} +tasks: + dtm: true + dsm: true + dhm: true + defaults: - _self_ # for hydra legacy reasons - - buffer: test.yaml - - filter: default_test.yaml + - interpolation: test.yaml + - dhm: test.yaml - io: test.yaml - tile_geometry: test.yaml # describes input features and classes - - dhm: test.yaml + - extract_stat: default.yaml # disable hydra logging - override hydra/hydra_logging: disabled diff --git a/environment.yml b/environment.yml index 3a12250..d545b76 100755 --- a/environment.yml +++ b/environment.yml @@ -6,11 +6,9 @@ dependencies: - python==3.10.* - pip - pytest - - cgal - gdal - laspy - numpy - - scipy - fiona - pyproj - pdal>=2.9.0 @@ -31,5 +29,4 @@ dependencies: - pip: - ign-pdal-tools - rasterstats - # - rasterio # moved here since pdal>=2.9.0 force modern libs to be used and make dependancy graph non resolved - + - rasterio # moved here since pdal>=2.9.0 force modern libs to be used and make dependancy graph non resolved \ No newline at end of file diff --git a/las_digital_models/add_buffer_one_tile.py b/las_digital_models/add_buffer_one_tile.py deleted file mode 100755 index 799d45d..0000000 --- a/las_digital_models/add_buffer_one_tile.py +++ /dev/null @@ -1,54 +0,0 @@ -"""Add a buffer around the queried tile from its neighbors -The script assumes that the neighbor tiles are located in the same folder as -the queried tile - -""" - -import logging -import os - -import hydra -from omegaconf import DictConfig -from pdaltools.las_add_buffer import create_las_with_buffer - -from las_digital_models.commons import commons - -log = commons.get_logger(__name__) - - -@hydra.main(config_path="../configs/", config_name="config.yaml", version_base="1.2") -def run_add_buffer_one_tile(config: DictConfig): - """ - The script assumes that the neighbor tiles are located in the same folder as - the queried tile - """ - - if config.io.forced_intermediate_ext is None: - input_file = os.path.join(config.io.input_dir, config.io.input_filename) - output_file = os.path.join(config.io.output_dir, config.io.input_filename) - else: - _, input_basename = os.path.split(config.io.input_filename) - tilename, _ = os.path.splitext(input_basename) - input_file = os.path.join(config.io.input_dir, f"{tilename}.{config.io.forced_intermediate_ext}") - output_file = os.path.join(config.io.output_dir, f"{tilename}.{config.io.forced_intermediate_ext}") - - os.makedirs(config.io.output_dir, exist_ok=True) - - create_las_with_buffer( - input_dir=config.io.input_dir, - tile_filename=input_file, - output_filename=output_file, - buffer_width=config.buffer.size, - spatial_ref=config.io.spatial_reference, - tile_width=config.tile_geometry.tile_width, - tile_coord_scale=config.tile_geometry.tile_coord_scale, - ) - - -def main(): - logging.basicConfig(level=logging.INFO) - run_add_buffer_one_tile() - - -if __name__ == "__main__": - main() diff --git a/las_digital_models/dhm_one_tile.py b/las_digital_models/dhm_one_tile.py index 889ecd2..d174451 100755 --- a/las_digital_models/dhm_one_tile.py +++ b/las_digital_models/dhm_one_tile.py @@ -12,26 +12,41 @@ log = commons.get_logger(__name__) -@hydra.main(config_path="../configs/", config_name="config.yaml", version_base="1.2") -def run_dhm_on_tile(config: DictConfig): - os.makedirs(config.io.output_dir, exist_ok=True) - tilename, _ = os.path.splitext(config.io.input_filename) +def run_dhm_on_tile( + input_las_filename: str, + input_dtm_dir: str, + input_dsm_dir: str, + output_dir: str, + pixel_size: float, + no_data_value: float, +): + os.makedirs(output_dir, exist_ok=True) + tilename, _ = os.path.splitext(input_las_filename) # for export - _size = commons.give_name_resolution_raster(config.tile_geometry.pixel_size) + _size = commons.give_name_resolution_raster(pixel_size) geotiff_filename = f"{tilename}{_size}.tif" - geotiff_dsm = os.path.join(config.dhm.input_dsm_dir, geotiff_filename) - geotiff_dtm = os.path.join(config.dhm.input_dtm_dir, geotiff_filename) - geotiff_output = os.path.join(config.io.output_dir, geotiff_filename) + geotiff_dsm = os.path.join(input_dsm_dir, geotiff_filename) + geotiff_dtm = os.path.join(input_dtm_dir, geotiff_filename) + geotiff_output = os.path.join(output_dir, geotiff_filename) # process - calculate_dhm(geotiff_dsm, geotiff_dtm, geotiff_output, no_data_value=config.tile_geometry.no_data_value) + calculate_dhm(geotiff_dsm, geotiff_dtm, geotiff_output, no_data_value=no_data_value) return -def main(): +@hydra.main(config_path="../configs/", config_name="config.yaml", version_base="1.2") +def main(config: DictConfig): logging.basicConfig(level=logging.INFO) - run_dhm_on_tile() + + run_dhm_on_tile( + input_las_filename=config.io.input_filename, + input_dtm_dir=config.dhm.input_dtm_dir, + input_dsm_dir=config.dhm.input_dsm_dir, + output_dir=config.io.output_dir, + pixel_size=config.tile_geometry.pixel_size, + no_data_value=config.tile_geometry.no_data_value, + ) if __name__ == "__main__": diff --git a/las_digital_models/ip_one_tile.py b/las_digital_models/ip_one_tile.py index 898a47b..e521cf4 100644 --- a/las_digital_models/ip_one_tile.py +++ b/las_digital_models/ip_one_tile.py @@ -1,4 +1,4 @@ -""" Main script for interpolation on a single tile +"""Main script for interpolation on a single tile Output files will be written to the target folder, tagged with the name of the interpolation method that was used. """ @@ -6,57 +6,100 @@ import logging import os import tempfile +from typing import List, Tuple import hydra from omegaconf import DictConfig +from pdaltools.las_info import parse_filename from las_digital_models.commons import commons -from las_digital_models.tasks.las_interpolation import interpolate_from_config +from las_digital_models.tasks.las_interpolation import interpolate from las_digital_models.tasks.postprocessing import mask_with_no_data_shapefile log = commons.get_logger(__name__) -@hydra.main(config_path="../configs/", config_name="config.yaml", version_base="1.2") -def run_ip_on_tile(config: DictConfig): - """Run interpolation on single tile using hydra config - config parameters are explained in the default.yaml files - """ - if config.io.input_dir is None: - input_dir = config.io.output_dir - else: - input_dir = config.io.input_dir - - os.makedirs(config.io.output_dir, exist_ok=True) - tilename, _ = os.path.splitext(config.io.input_filename) - - # input file (already filtered and potentially with a buffer) - if config.io.forced_intermediate_ext is None: - input_file = os.path.join(input_dir, config.io.input_filename) - else: - input_file = os.path.join(input_dir, f"{tilename}.{config.io.forced_intermediate_ext}") +def run_ip_on_tile( + input_dir: str, + input_filename: str, + output_dir: str, + origin: Tuple[int, int], + pixel_size: float, + tile_width: int, + spatial_reference: str, + no_data_value: float, + no_data_mask_shapefile: str, + filter_dimension: str, + filter_keep_values: List[int], +): + """Run interpolation on single tile with geometry masking in case a mask is provided""" - # for export - _size = commons.give_name_resolution_raster(config.tile_geometry.pixel_size) + # Generate output filenames + tilename, _ = os.path.splitext(input_filename) + _size = commons.give_name_resolution_raster(pixel_size) geotiff_stem = f"{tilename}{_size}" geotiff_filename = f"{geotiff_stem}.tif" - geotiff_path = os.path.join(config.io.output_dir, geotiff_filename) + geotiff_path = os.path.join(output_dir, geotiff_filename) - if config.io.no_data_mask_shapefile: + input_path = os.path.join(input_dir, input_filename) + os.makedirs(output_dir, exist_ok=True) + + if no_data_mask_shapefile: with tempfile.NamedTemporaryFile(suffix=".tif", prefix=f"{geotiff_stem}_raw") as tmp_geotiff: # process interpolation - interpolate_from_config(input_file, tmp_geotiff.name, config) - mask_with_no_data_shapefile( - config.io.no_data_mask_shapefile, tmp_geotiff.name, geotiff_path, config.tile_geometry.no_data_value + interpolate( + input_path, + tmp_geotiff.name, + origin, + pixel_size, + tile_width, + spatial_reference, + no_data_value, + filter_dimension, + filter_keep_values, ) + mask_with_no_data_shapefile(no_data_mask_shapefile, tmp_geotiff.name, geotiff_path, no_data_value) else: - interpolate_from_config(input_file, geotiff_path, config) + interpolate( + input_path, + geotiff_path, + origin, + pixel_size, + tile_width, + spatial_reference, + no_data_value, + filter_dimension, + filter_keep_values, + ) -def main(): +@hydra.main(config_path="../configs/", config_name="config.yaml", version_base="1.2") +def main(config: DictConfig): logging.basicConfig(level=logging.INFO) - run_ip_on_tile() + + # Use filename to get the tile coordinates. + # Coordinates are needed to define neighboring tiles, in order to create the buffered tile. + input_path = os.path.join(config.io.input_dir, config.io.input_filename) + _, coordX, coordY, _ = parse_filename(input_path) + origin = [ + float(coordX) * config.tile_geometry.tile_coord_scale, + float(coordY) * config.tile_geometry.tile_coord_scale, + ] + + run_ip_on_tile( + input_dir=config.io.input_dir, + input_filename=config.io.input_filename, + origin=origin, + output_dir=config.io.output_dir, + pixel_size=config.tile_geometry.pixel_size, + tile_width=config.tile_geometry.tile_width, + spatial_reference=config.io.spatial_reference, + no_data_value=config.tile_geometry.no_data_value, + no_data_mask_shapefile=config.io.no_data_mask_shapefile, + filter_dimension=config.interpolation.custom.filter.dimension, + filter_keep_values=config.interpolation.custom.filter.keep_values, + ) if __name__ == "__main__": diff --git a/las_digital_models/main.py b/las_digital_models/main.py new file mode 100644 index 0000000..dc5218e --- /dev/null +++ b/las_digital_models/main.py @@ -0,0 +1,140 @@ +"""Main entry point for digital models generation""" + +import logging as log +import os +import tempfile +from pathlib import Path + +import hydra +from omegaconf import DictConfig +from pdaltools.las_add_buffer import create_las_with_buffer +from pdaltools.las_info import get_tile_origin_using_header_info + +from las_digital_models.dhm_one_tile import run_dhm_on_tile +from las_digital_models.ip_one_tile import run_ip_on_tile + + +@hydra.main(config_path="../configs/", config_name="config.yaml", version_base="1.2") +def main_las_digital_models(config: DictConfig): + """Main function to generate digital models from a las file. + It can compute: + - DSM: digital surface model + - DTM: digital terrain model + - DHM: digital height model (DSM - DTM) + + If a buffer size is set in config, the las file is buffered using its neighbors + before generating the digital models to prevent border effects. + Neighbors search is computed based on file names: + las files are expected to be formatted as: {prefix1}_{prefix2}_{XXXX}_{YYYY}_{suffix} + with XXXX and YYYY their coordinates with conditions: + - XXXX (and YYYY) * config.tile_geometry.tile_coord_scale are the coordinates of the upper left + corner in meters + - XXXX (and YYYY) * config.tile_geometry.tile_coord_scale + (- for YYYY) config.tile_geometry.tile_width + are the coordinates of the lower right corner in meters + + Args: + config (DictConfig): hydra config for the project + """ + log.basicConfig(level=log.INFO, format="%(message)s") + + initial_las_filename = config.io.input_filename + in_dir = Path(config.io.input_dir) + out_dir = Path(config.io.output_dir) + + # Check input/output files and folders + if initial_las_filename is None or in_dir is None or out_dir is None: + raise RuntimeError( + """In input you have to give a las, an input directory and an output directory. + For more info run the same command by adding --help""" + ) + + os.makedirs(out_dir, exist_ok=True) + + with ( + tempfile.TemporaryDirectory(prefix="tmp_buffer", dir=".") as tmpdir_buffer, + tempfile.TemporaryDirectory(prefix="tmp_dtm", dir=".") as tmpdir_dtm, + tempfile.TemporaryDirectory(prefix="tmp_dsm", dir=".") as tmpdir_dsm, + ): + # Get pointcloud origin from the las file metadata + tile_origin = get_tile_origin_using_header_info( + in_dir / initial_las_filename, tile_width=config.tile_geometry.tile_width + ) + # Buffer + if config.buffer.size: + log.info(f"Create buffered las file with buffer = {config.buffer.size}") + if config.buffer.output_subdir: + buffer_output_dir = Path(out_dir) / config.buffer.output_subdir + else: + buffer_output_dir = Path(tmpdir_buffer) + buffer_output_dir.mkdir(parents=True, exist_ok=True) + + epsg = config.io.spatial_reference + create_las_with_buffer( + input_dir=str(in_dir), + tile_filename=str(in_dir / initial_las_filename), + output_filename=str(buffer_output_dir / initial_las_filename), + buffer_width=config.buffer.size, + spatial_ref=f"EPSG:{epsg}" if str(epsg).isdigit() else epsg, + tile_width=config.tile_geometry.tile_width, + tile_coord_scale=config.tile_geometry.tile_coord_scale, + ) + else: + log.info("Skip las buffer creation") + buffer_output_dir = in_dir + + # Compute DTM + if config.tasks.dtm or config.tasks.dhm: + log.info("Create DTM") + dtm_output_dir = ( + os.path.join(out_dir, config.interpolation.dtm.output_subfolder) if config.tasks.dtm else tmpdir_dtm + ) + run_ip_on_tile( + input_dir=buffer_output_dir, + input_filename=initial_las_filename, + origin=tile_origin, + output_dir=dtm_output_dir, + pixel_size=config.tile_geometry.pixel_size, + tile_width=config.tile_geometry.tile_width, + spatial_reference=config.io.spatial_reference, + no_data_value=config.tile_geometry.no_data_value, + no_data_mask_shapefile=config.io.no_data_mask_shapefile, + filter_dimension=config.interpolation.dtm.filter.dimension, + filter_keep_values=config.interpolation.dtm.filter.keep_values, + ) + + # Compute DSM + if config.tasks.dsm or config.tasks.dhm: + log.info("Create DSM") + dsm_output_dir = ( + os.path.join(out_dir, config.interpolation.dsm.output_subfolder) if config.tasks.dsm else tmpdir_dsm + ) + run_ip_on_tile( + input_dir=buffer_output_dir, + input_filename=initial_las_filename, + origin=tile_origin, + output_dir=dsm_output_dir, + pixel_size=config.tile_geometry.pixel_size, + tile_width=config.tile_geometry.tile_width, + spatial_reference=config.io.spatial_reference, + no_data_value=config.tile_geometry.no_data_value, + no_data_mask_shapefile=config.io.no_data_mask_shapefile, + filter_dimension=config.interpolation.dsm.filter.dimension, + filter_keep_values=config.interpolation.dsm.filter.keep_values, + ) + + # Compute DHM + if config.tasks.dhm: + log.info("Create DHM") + dhm_output_dir = os.path.join(out_dir, config.dhm.output_subfolder) + run_dhm_on_tile( + input_las_filename=initial_las_filename, + input_dtm_dir=dtm_output_dir, + input_dsm_dir=dsm_output_dir, + output_dir=dhm_output_dir, + pixel_size=config.tile_geometry.pixel_size, + no_data_value=config.tile_geometry.no_data_value, + ) + + +if __name__ == "__main__": + main_las_digital_models() diff --git a/las_digital_models/tasks/las_interpolation.py b/las_digital_models/tasks/las_interpolation.py index 9ef070c..5e31966 100755 --- a/las_digital_models/tasks/las_interpolation.py +++ b/las_digital_models/tasks/las_interpolation.py @@ -1,60 +1,20 @@ -from typing import List +from typing import List, Tuple import pdal from osgeo import gdal -from pdaltools.las_info import parse_filename from las_digital_models.commons import commons gdal.UseExceptions() -def interpolate_from_config(input_file: str, output_raster: str, config: dict): - """API using a config dictionary for the `interpolate` method defined in this file - Generate a Z (height) raster file from a LAS point cloud file by interpolating the Z value for each pixel center. - - - Args: - input_file (str): path to the las/laz file to interpolate - output_file (str): path to the output raster - config (dict): ProduitDeriveLidar config dictionary containing - { - "tile_geometry": { - "tile_coord_scale": #int, scale of the tiles coordinates in the las filename - "tile_width": #int, width of the tile in meters (used to infer the lower-left corner) - "pixel_size": #float, pixel size of the output raster in meters (pixels are supposed to be squares) - "no_data_value": #int, no data value for the output raster - }, - "io": { - "spatial_reference": #str, spatial reference to use when reading las file - }, - "filter": { - "dimension": #str, dimension alogn which to filter - "keep_values": #list of ints, values of the gilter dimension for the points to use in the interpolation - } - } - - """ - interpolate( - input_file, - output_raster, - config["tile_geometry"]["pixel_size"], - config["tile_geometry"]["tile_width"], - config["tile_geometry"]["tile_coord_scale"], - config["io"]["spatial_reference"], - config["tile_geometry"]["no_data_value"], - config["filter"]["dimension"], - config["filter"]["keep_values"], - ) - - @commons.eval_time_with_pid def interpolate( input_file: str, output_file: str, + tile_origin: Tuple[int, int], pixel_size: float, tile_width: int, - tile_coord_scale: int, spatial_ref: str, no_data_value: int, filter_dimension: str, @@ -75,9 +35,9 @@ def interpolate( Args: input_file (str): path to the las/laz file to interpolate output_file (str): path to the output raster + tile_origin (Tuple[int, int]): coordinates of the upper left corner of the tile (in meters) pixel_size (float): pixel size of the output raster in meters (pixels are supposed to be squares) tile_width (int): width of the tile in meters (used to infer the lower-left corner) - tile_coord_scale (int): scale of the tiles coordinates in the las filename spatial_ref (str): spatial reference to use when reading las file no_data_value (int): no data value for the output raster filter_dimension (str): Name of the dimension along which to filter input points @@ -85,10 +45,6 @@ def interpolate( filter_values (List[int]): Values to keep for input points along filter_dimension """ - _, coordX, coordY, _ = parse_filename(input_file) - - # Compute origin/number of pixels - origin = [float(coordX) * tile_coord_scale, float(coordY) * tile_coord_scale] nb_pixels = [int(tile_width / pixel_size), int(tile_width / pixel_size)] # Read with pdal @@ -101,8 +57,8 @@ def interpolate( pipeline |= pdal.Filter.faceraster( resolution=str(pixel_size), - origin_x=str(origin[0] - pixel_size / 2), # lower left corner - origin_y=str(origin[1] + pixel_size / 2 - tile_width), # lower left corner + origin_x=str(tile_origin[0] - pixel_size / 2), # lower left corner + origin_y=str(tile_origin[1] + pixel_size / 2 - tile_width), # lower left corner width=str(nb_pixels[0]), height=str(nb_pixels[1]), ) diff --git a/run.sh b/run.sh index d9894a8..c846556 100755 --- a/run.sh +++ b/run.sh @@ -48,80 +48,13 @@ FILENAMES=($(cd ${INPUT} && find . -maxdepth 1 -type f \( -iname \*.las -o -ina echo "GENERATE DSM/DTM/DHM ON ${#FILENAMES[@]} FILES" echo "" -# Step 0: Create las with buffer from its neighbors tiles -# /!\ rasters generated from these las tiles will need to be cropped to match the input las area -BUFFERED_DIR=${OUTPUT}/las_with_buffer - -echo "------------------" -echo "add buffer" -echo "------------------" parallel --jobs ${PARALLEL_JOBS} \ - python -m las_digital_models.add_buffer_one_tile \ + python -m las_digital_models.main \ --config-name=${CONFIG_NAME} \ io.input_dir=${INPUT} \ io.input_filename={} \ - io.output_dir=${BUFFERED_DIR} \ - ::: "${FILENAMES[@]}" - - - -echo "------------------" -echo "RUN DTM GENERATION" -echo "------------------" - -# Output filenames for each step -DTM_DIR=${OUTPUT}/DTM - - -# Step 1 ; create DTM -echo "------------------" -echo "interpolation" -echo "------------------" - -parallel --jobs ${PARALLEL_JOBS} \ - python -m las_digital_models.ip_one_tile \ - --config-name=${CONFIG_NAME} \ - io.input_dir=${BUFFERED_DIR} \ - io.input_filename={} \ - filter=dtm \ - io.output_dir=${DTM_DIR} \ + io.output_dir=${OUTPUT} \ tile_geometry.pixel_size=${PIXEL_SIZE} \ io.no_data_mask_shapefile=${SHAPEFILE} \ ::: "${FILENAMES[@]}" - -echo "------------------" -echo "RUN DSM GENERATION" -echo "------------------" - -DSM_DIR=${OUTPUT}/DSM - -# Step 2; create DSM -parallel --jobs ${PARALLEL_JOBS} \ - python -m las_digital_models.ip_one_tile \ - --config-name=${CONFIG_NAME} \ - io.input_dir=${BUFFERED_DIR} \ - io.input_filename={} \ - filter=dsm \ - io.output_dir=${DSM_DIR} \ - tile_geometry.pixel_size=${PIXEL_SIZE} \ - io.no_data_mask_shapefile=${SHAPEFILE} \ - ::: "${FILENAMES[@]}" - -echo "------------------" -echo "RUN DHM GENERATION" -echo "------------------" - -# Step 3 ; create DHM with each of the methods -# Output filenames for each step -DHM_DIR=${OUTPUT}/DHM - -parallel --jobs ${PARALLEL_JOBS} \ - python -m las_digital_models.dhm_one_tile \ - --config-name=${CONFIG_NAME} \ - dhm.input_dsm_dir=${DSM_DIR} \ - dhm.input_dtm_dir=${DTM_DIR} \ - io.input_filename={} \ - io.output_dir=${DHM_DIR} \ - tile_geometry.pixel_size=${PIXEL_SIZE} \ - ::: "${FILENAMES[@]}" diff --git a/test/tasks/test_las_interpolation.py b/test/tasks/test_las_interpolation.py index 96743dd..27c9b57 100644 --- a/test/tasks/test_las_interpolation.py +++ b/test/tasks/test_las_interpolation.py @@ -18,9 +18,10 @@ COORD_X = 77055 COORD_Y = 627760 -EXPECTED_XMIN = COORD_X * TILE_COORD_SCALE - PIXEL_SIZE / 2 -EXPECTED_XMAX = COORD_Y * TILE_COORD_SCALE + PIXEL_SIZE / 2 -EXPECTED_RASTER_BOUNDS = (EXPECTED_XMIN, EXPECTED_XMAX - TILE_WIDTH), (EXPECTED_XMIN + TILE_WIDTH, EXPECTED_XMAX) +ORIGIN = [COORD_X * TILE_COORD_SCALE, COORD_Y * TILE_COORD_SCALE] +EXPECTED_XMIN = ORIGIN[0] - PIXEL_SIZE / 2 +EXPECTED_YMAX = ORIGIN[1] + PIXEL_SIZE / 2 +EXPECTED_RASTER_BOUNDS = (EXPECTED_XMIN, EXPECTED_YMAX - TILE_WIDTH), (EXPECTED_XMIN + TILE_WIDTH, EXPECTED_YMAX) def setup_module(): @@ -62,9 +63,9 @@ def test_interpolate(filter_dimension, filter_values, output_file, ground_truth_ interpolate( INPUT_FILE, output_file, + tile_origin=ORIGIN, pixel_size=PIXEL_SIZE, tile_width=TILE_WIDTH, - tile_coord_scale=TILE_COORD_SCALE, spatial_ref="EPSG:2154", no_data_value=-9999, filter_dimension=filter_dimension, diff --git a/test/test_add_buffer_one_tile.py b/test/test_add_buffer_one_tile.py deleted file mode 100644 index 7ba669e..0000000 --- a/test/test_add_buffer_one_tile.py +++ /dev/null @@ -1,55 +0,0 @@ -import logging -import os -import shutil -import test.utils.point_cloud_utils as pcu - -from hydra import compose, initialize - -from las_digital_models import add_buffer_one_tile - -TEST_PATH = os.path.dirname(__file__) -TMP_PATH = os.path.join(TEST_PATH, "tmp/buffer") -DATA_PATH = os.path.join(TEST_PATH, "data") - - -def setup_module(module): - try: - shutil.rmtree(TMP_PATH) - - except FileNotFoundError: - pass - os.makedirs(TMP_PATH, exist_ok=True) - - -def test_add_buffer_one_tile(): - # DATA_PATH contains a .laz file, but check that the io.forced_intermediate_ext parameter forces to look for a - # .laz file - input_filename = "test_data_77055_627760_LA93_IGN69.las" - output_filename = "test_data_77055_627760_LA93_IGN69.laz" - output_file = os.path.join(TMP_PATH, output_filename) - - with initialize(version_base="1.2", config_path="../configs"): - # config is relative to a module - cfg = compose( - config_name="config", - overrides=[ - "io=test", - "tile_geometry=test", - f"io.input_dir={DATA_PATH}", - f"io.input_filename={input_filename}", - "io.forced_intermediate_ext=laz", - f"io.output_dir={TMP_PATH}", - "buffer=test", - ], - ) - - add_buffer_one_tile.run_add_buffer_one_tile(cfg) - assert os.path.isfile(output_file) - logging.info(pcu.get_nb_points(output_file)) - assert pcu.get_nb_points(output_file) == 103359 - assert pcu.get_classification_values(output_file) == {1, 2, 3, 4, 5, 6, 64} - - -if __name__ == "__main__": - logging.basicConfig(level=logging.INFO) - test_add_buffer_one_tile() diff --git a/test/test_dhm_one_tile.py b/test/test_dhm_one_tile.py index fb3a439..3c454c7 100644 --- a/test/test_dhm_one_tile.py +++ b/test/test_dhm_one_tile.py @@ -1,4 +1,3 @@ -import logging import os import shutil import test.utils.raster_utils as ru @@ -33,7 +32,7 @@ def setup_module(module): os.mkdir(tmp_path) -def test_mnh_one_tile(): +def test_dhm_one_tile(): with initialize(version_base="1.2", config_path="../configs"): # config is relative to a module cfg = compose( @@ -41,18 +40,13 @@ def test_mnh_one_tile(): overrides=[ "io=test", "tile_geometry=test", - f"io.output_dir={output_dir}", "dhm=test", + f"io.output_dir={output_dir}", ], ) - dhm_one_tile.run_dhm_on_tile(cfg) + dhm_one_tile.main(cfg) assert os.path.isfile(expected_output_file) raster_bounds = ru.get_tif_extent(expected_output_file) assert ru.allclose_mm(raster_bounds, expected_raster_bounds) - - -if __name__ == "__main__": - logging.basicConfig(level=logging.DEBUG) - test_mnh_one_tile() diff --git a/test/test_ip_one_tile.py b/test/test_ip_one_tile.py index 9aedf2c..bd74081 100644 --- a/test/test_ip_one_tile.py +++ b/test/test_ip_one_tile.py @@ -15,12 +15,12 @@ PIXEL_SIZE = 0.5 TEST_PATH = Path(__file__).resolve().parent -TMP_PATH = TEST_PATH / "tmp" +TMP_PATH = TEST_PATH / "tmp" / "ip_one_tile" GROUND_TRUTH_FOLDER = TEST_PATH / "data" / "interpolation" EXPECTED_XMIN = COORD_X * TILE_COORD_SCALE - PIXEL_SIZE / 2 -EXPECTED_XMAX = COORD_Y * TILE_COORD_SCALE + PIXEL_SIZE / 2 -EXPECTED_RASTER_BOUNDS = (EXPECTED_XMIN, EXPECTED_XMAX - TILE_WIDTH), (EXPECTED_XMIN + TILE_WIDTH, EXPECTED_XMAX) +EXPECTED_YMAX = COORD_Y * TILE_COORD_SCALE + PIXEL_SIZE / 2 +EXPECTED_RASTER_BOUNDS = (EXPECTED_XMIN, EXPECTED_YMAX - TILE_WIDTH), (EXPECTED_XMIN + TILE_WIDTH, EXPECTED_YMAX) SHAPEFILE = TEST_PATH / "data" / "mask_shapefile" / "test_multipolygon_shapefile.shp" EXPECTED_OUTPUT_USING_SHAPEFILE = GROUND_TRUTH_FOLDER / "test_data_77055_627760_LA93_IGN69_50CM_no_data.tif" @@ -43,8 +43,8 @@ def get_expected_output_file(base_dir=None): return expected_output_file -def test_ip_one_tile(): - output_dir = os.path.join(TMP_PATH, "test_ip_one_tile") +def test_main(): + output_dir = os.path.join(TMP_PATH, "test_main") os.makedirs(output_dir, exist_ok=True) with initialize(version_base="1.2", config_path="../configs"): # config is relative to a module @@ -54,15 +54,15 @@ def test_ip_one_tile(): "io=test", f"io.output_dir={output_dir}", "tile_geometry=test", - "filter.dimension=''", - "filter.keep_values=[]", + "interpolation.custom.filter.dimension=''", + "interpolation.custom.filter.keep_values=[]", ], ) output_file = get_expected_output_file(base_dir=output_dir) logging.debug(output_file) logging.debug(f"Pixel size: {cfg.tile_geometry.pixel_size}") - ip_one_tile.run_ip_on_tile(cfg) + ip_one_tile.main(cfg) assert os.path.isfile(output_file) raster_bounds = ru.get_tif_extent(output_file) @@ -71,8 +71,8 @@ def test_ip_one_tile(): assert ru.tif_values_all_close(output_file, os.path.join(GROUND_TRUTH_FOLDER, os.path.basename(output_file))) -def test_ip_with_no_data_mask(): - output_dir = os.path.join(TMP_PATH, "test_ip_with_no_data_mask") +def test_main_with_no_data_mask(): + output_dir = os.path.join(TMP_PATH, "test_main_with_no_data_mask") os.makedirs(output_dir, exist_ok=True) with initialize(version_base="1.2", config_path="../configs"): # config is relative to a module @@ -83,52 +83,17 @@ def test_ip_with_no_data_mask(): "tile_geometry=test", f"io.no_data_mask_shapefile={SHAPEFILE}", f"io.output_dir={output_dir}", - "filter.keep_values=[]", + "interpolation.custom.filter.keep_values=[]", ], ) output_file = get_expected_output_file(base_dir=output_dir) logging.debug(f"Write to {output_file}") - ip_one_tile.run_ip_on_tile(cfg) + ip_one_tile.main(cfg) assert os.path.isfile(output_file) raster_bounds = ru.get_tif_extent(output_file) assert ru.allclose_mm(raster_bounds, EXPECTED_RASTER_BOUNDS) assert ru.tif_values_all_close(output_file, EXPECTED_OUTPUT_USING_SHAPEFILE) - - -def test_ip_dtm_classes(): - output_dir = os.path.join(TMP_PATH, "test_ip_dtm_classes") - os.makedirs(output_dir, exist_ok=True) - expected_output = os.path.join( - TEST_PATH, "data", "interpolation", "test_data_77055_627760_LA93_IGN69_50CM_dtm_classes.tif" - ) - with initialize(version_base="1.2", config_path="../configs"): - # config is relative to a module - cfg = compose( - config_name="config", - overrides=[ - "io=test", - "tile_geometry=test", - "filter=dtm", - f"io.output_dir={output_dir}", - ], - ) - - output_file = get_expected_output_file(base_dir=output_dir) - - ip_one_tile.run_ip_on_tile(cfg) - assert os.path.isfile(output_file) - - raster_bounds = ru.get_tif_extent(output_file) - assert ru.allclose_mm(raster_bounds, EXPECTED_RASTER_BOUNDS) - - assert ru.tif_values_all_close(output_file, expected_output) - - -if __name__ == "__main__": - logging.basicConfig(level=logging.DEBUG) - test_ip_one_tile() - test_ip_with_no_data_mask() diff --git a/test/test_main.py b/test/test_main.py new file mode 100644 index 0000000..fbf30aa --- /dev/null +++ b/test/test_main.py @@ -0,0 +1,131 @@ +import os +import shutil +import test.utils.raster_utils as ru +from pathlib import Path + +import laspy +import numpy as np +from hydra import compose, initialize + +from las_digital_models.main import main_las_digital_models + +COORD_X = 77055 +COORD_Y = 627760 +TILE_COORD_SCALE = 10 +TILE_WIDTH = 50 +PIXEL_SIZE = 0.5 +TEST_PATH = Path(__file__).resolve().parent +TMP_PATH = TEST_PATH / "tmp" / "main" +DATA_PATH = TEST_PATH / "data" + +INPUT_FILENAME = f"test_data_{COORD_X}_{COORD_Y}_LA93_IGN69.laz" +INPUT_TILENAME = os.path.splitext(INPUT_FILENAME)[0] +OUTPUT_TIF_NAME = f"{INPUT_TILENAME}_50CM.tif" + +expected_xmin = COORD_X * TILE_COORD_SCALE - PIXEL_SIZE / 2 +expected_ymax = COORD_Y * TILE_COORD_SCALE + PIXEL_SIZE / 2 +EXPECTED_RASTER_BOUNDS = (expected_xmin, expected_ymax - TILE_WIDTH), (expected_xmin + TILE_WIDTH, expected_ymax) + + +def setup_module(): + try: + shutil.rmtree(TMP_PATH) + + except FileNotFoundError: + pass + os.makedirs(TMP_PATH) + + +def get_2d_bounding_box(path): + """Get bbox for a las file (x, y only)""" + with laspy.open(path) as f: + mins = f.header.mins + maxs = f.header.maxs + + return mins[:2], maxs[:2] + + +def test_main_intermediate_files(): + in_mins, in_maxs = get_2d_bounding_box(DATA_PATH / INPUT_FILENAME) + buffer_size = 10 + output_dir = TMP_PATH / "main_intermediate_files" + output_buffer_dir = output_dir / "buffer" + output_dtm_dir = output_dir / "DTM" + output_dsm_dir = output_dir / "DSM" + output_dhm_dir = output_dir / "DHM" + + with initialize(version_base="1.2", config_path="../configs"): + # config is relative to a module + cfg = compose( + config_name="test", + overrides=[ + f"io.input_filename={INPUT_FILENAME}", + f"io.input_dir={DATA_PATH}", + f"io.output_dir={output_dir}", + f"tile_geometry.tile_coord_scale={TILE_COORD_SCALE}", + f"tile_geometry.tile_width={TILE_WIDTH}", + f"buffer.size={buffer_size}", + "buffer.output_subdir='buffer'", + ], + ) + main_las_digital_models(cfg) + + # Check buffer files are correct + output_buffer_path = output_buffer_dir / INPUT_FILENAME + assert os.path.isfile(output_buffer_path) + out_mins, out_maxs = get_2d_bounding_box(output_buffer_path) + assert np.all(out_mins == in_mins - buffer_size) + assert np.all(out_maxs[0] == in_maxs[0] + buffer_size) + assert out_maxs[1] == in_maxs[1] # neighbor file does not exist + + # Check output tif files are correct + output_dtm_path = output_dtm_dir / OUTPUT_TIF_NAME + assert os.path.isfile(output_dtm_path) + dtm_bounds = ru.get_tif_extent(output_dtm_path) + + assert dtm_bounds == EXPECTED_RASTER_BOUNDS + + output_dsm_path = output_dsm_dir / OUTPUT_TIF_NAME + assert os.path.isfile(output_dsm_path) + dsm_bounds = ru.get_tif_extent(output_dsm_path) + + assert dsm_bounds == EXPECTED_RASTER_BOUNDS + + output_dhm_path = output_dhm_dir / OUTPUT_TIF_NAME + assert os.path.isfile(output_dhm_path) + dhm_bounds = ru.get_tif_extent(output_dhm_path) + + assert dhm_bounds == EXPECTED_RASTER_BOUNDS + + +def test_main_without_intermediate_files(): + """Compute only dhm, and check buffer / dsm / dtm files are only temporary files""" + buffer_size = 10 + output_dir = TMP_PATH / "main_without_intermediate_files" + output_buffer_dir = output_dir / "buffer" + output_dtm_dir = output_dir / "DTM" + output_dsm_dir = output_dir / "DSM" + output_dhm_dir = output_dir / "DHM" + with initialize(version_base="1.2", config_path="../configs"): + # config is relative to a module + cfg = compose( + config_name="test", + overrides=[ + f"io.input_filename={INPUT_FILENAME}", + f"io.input_dir={DATA_PATH}", + f"io.output_dir={output_dir}", + f"tile_geometry.tile_coord_scale={TILE_COORD_SCALE}", + f"tile_geometry.tile_width={TILE_WIDTH}", + f"buffer.size={buffer_size}", + "tasks.dtm=false", + "tasks.dsm=false", + "tasks.dhm=true", + ], + ) + main_las_digital_models(cfg) + + # Check outputs are correct + assert not os.path.exists(output_buffer_dir) + assert not os.path.exists(output_dtm_dir) + assert not os.path.exists(output_dsm_dir) + assert os.path.exists(output_dhm_dir / OUTPUT_TIF_NAME) diff --git a/test/test_run_script.py b/test/test_run_script.py index 7073ac5..57ff608 100644 --- a/test/test_run_script.py +++ b/test/test_run_script.py @@ -8,9 +8,9 @@ from las_digital_models.commons import commons test_path = os.path.dirname(os.path.abspath(__file__)) -tmp_path = os.path.join(test_path, "tmp") +tmp_path = os.path.join(test_path, "tmp", "run_script") input_dir = os.path.join(test_path, "data") -output_dir = os.path.join(tmp_path, "output_run_script") +output_dir = tmp_path file_ext = "laz" pixel_size = 0.5 @@ -27,7 +27,7 @@ def setup_module(module): except FileNotFoundError: pass - os.mkdir(tmp_path) + os.makedirs(tmp_path) @pytest.mark.functional_test