From ae9016477f3a8f81e436b85694136c0f0328ac63 Mon Sep 17 00:00:00 2001 From: ShashankBice Date: Fri, 21 Feb 2025 17:30:13 +0000 Subject: [PATCH 01/11] add functionality to mosaic --- src/grid_pc/dsm_functions.py | 105 +++++++++++++++++++++++++++++++++++ src/grid_pc/pdal_pipeline.py | 55 +++++++++++++----- 2 files changed, 147 insertions(+), 13 deletions(-) diff --git a/src/grid_pc/dsm_functions.py b/src/grid_pc/dsm_functions.py index 86b5715..b9a27fb 100644 --- a/src/grid_pc/dsm_functions.py +++ b/src/grid_pc/dsm_functions.py @@ -6,6 +6,9 @@ from shapely.geometry import Polygon import geopandas as gpd import requests +import subprocess +import ast +from shutil import which def return_readers(input_aoi, src_crs, @@ -198,3 +201,105 @@ def create_dem_stage(dem_filename='dem_output.tif', pointcloud_resolution=1., return [dem_stage] +def dem_mosaic(img_list,outfn,tr=None,tsrs=None,stats=None,tile_size=None,extent=None): + """ + From https://github.com/uw-cryo/skysat_stereo/blob/master/skysat_stereo/asp_utils.py + mosaic input image list using ASP's dem_mosaic program. + See dem_mosaic documentation here: https://stereopipeline.readthedocs.io/en/latest/tools/dem_mosaic.html + Parameters + ---------- + img_list: list + List of input images to be mosaiced + outfn: str + Path to output mosaicked image + tr: float/int + target resolution of output mosaic + t_srs: str + target projection of output mosaic (default: EPSG:4326) + stats: str + metric to use for mosaicing + tile_size: int + tile size for distributed mosaicing (if less on memory) + Returns + ---------- + out: str + dem_mosaic log + """ + + dem_mosaic_opt = [] + + if stats: + dem_mosaic_opt.extend(['--{}'.format(stats)]) + if (tr is not None) & (ast.literal_eval(tr) is not None): + dem_mosaic_opt.extend(['--tr', str(tr)]) + if tsrs: + dem_mosaic_opt.extend(['--t_srs', tsrs]) + if extent: + xmin,ymin,xmax,ymax = extent.split(' ') + dem_mosaic_opt.extend(['--t_projwin', xmin,ymin,xmax,ymax]) + dem_mosaic_args = img_list + if tile_size: + # will first perform tile-wise vertical mosaicing + # then blend the result + dem_mosaic_opt.extend(['--tile-size',str(tile_size)]) + temp_fol = os.path.splitext(outfn)[0]+'_temp' + dem_mosaic_opt.extend(['-o',os.path.join(temp_fol,'run')]) + out_tile_op = run_cmd('dem_mosaic',dem_mosaic_args+dem_mosaic_opt) + # query all tiles and then do simple mosaic + #print(os.path.join(temp_fol,'run-*.tif')) + mos_tile_list = sorted(glob.glob(os.path.join(temp_fol,'run-*.tif'))) + print(f"Found {len(mos_tile_list)}") + # now perform simple mosaic + dem_mos2_opt = [] + dem_mos2_opt.extend(['-o',outfn]) + dem_mos2_args = mos_tile_list + out_fn_mos = run_cmd('dem_mosaic',dem_mos2_args+dem_mos2_opt) + out = out_tile_op+out_fn_mos + print("Deleting tile directory") + shutil.rmtree(temp_fol) + + else: + # process all at once, no tiling + dem_mosaic_opt.extend(['-o',outfn]) + out = run_cmd('dem_mosaic',dem_mosaic_args+dem_mosaic_opt) + return out + +def run_cmd(bin, args, **kw): + """ + From https://github.com/uw-cryo/skysat_stereo/blob/master/skysat_stereo/asp_utils.py + wrapper around subprocess function to excute bash commands + Parameters + ---------- + bin: str + command to be excuted (e.g., stereo or gdalwarp) + args: list + arguments to the command as a list + Retuns + ---------- + out: str + log (stdout) as str if the command executed, error message if the command failed + """ + + #from dshean/vmap.py + + binpath = which(bin) + #if binpath is None: + #msg = ("Unable to find executable %s\n" + #"Install ASP and ensure it is in your PATH env variable\n" + #"https://ti.arc.nasa.gov/tech/asr/intelligent-robotics/ngt/stereo/" % bin) + #sys.exit(msg) + #binpath = os.path.join('/opt/StereoPipeline/bin/',bin) + call = [binpath,] + if args is not None: + call.extend(args) + #print(call) + try: + out = subprocess.run(call,check=True,capture_output=True,encoding='UTF-8').stdout + except: + out = "the command {} failed to run, see corresponding asp log".format(call) + return out + + + + + diff --git a/src/grid_pc/pdal_pipeline.py b/src/grid_pc/pdal_pipeline.py index 4650b7b..a94323a 100644 --- a/src/grid_pc/pdal_pipeline.py +++ b/src/grid_pc/pdal_pipeline.py @@ -36,8 +36,9 @@ def create_dsm(extent_geojson: str, # processing_extent.geojson source_wkt: str | None, # SRS_CRS.wkt target_wkt: str, # UTM_13N_WGS84_G2139_3D.wkt - output_path: str, # /tmp/dem/ - ) -> None: + output_prefix: str, # CO_ALS_proc/CO_3DEP_ALS + mosaic: bool | True, + cleanup: bool | True) -> None: """ Create a Digital Surface Model (DSM) from a given extent and point cloud data. This function divides a region of interest into tiles and generates a DSM geotiff from EPT point clouds for each tile using PDAL @@ -50,18 +51,20 @@ def create_dsm(extent_geojson: str, # processing_extent.geojson Path to the WKT file defining the source coordinate reference system (CRS). If None, the CRS from the point cloud file is used. target_wkt : str Path to the WKT file defining the target coordinate reference system (CRS). - output_path : str - Directory path where the output DSM files will be saved. - + output_prefix : str + prefix with directory name and filename prefix for the project (e.g., CO_ALS_proc/CO_3DEP_ALS) + mosaic : bool + Mosaic the output tiles using a weighted average algorithm + cleanup: bool + If true, remove the intermediate tif files for the output tiles Returns ------- None Notes ----- - After running this function, create a final DEM with the following commands: - 1. dem_moaic -o merged_dsm.tif dem*.tif - 2. gdalwarp -s_srs SRS_CRS.wkt -t_srs UTM_13N_WGS84_G2139_3D.wkt -r cubic -tr 1.0 1.0 merged_dsm.tif merged_dsm_reprojected_UTM_13N_WGS84_G2139.tif + After running this function, reproject the final DEM with the following commands: + 1. gdalwarp -s_srs SRS_CRS.wkt -t_srs UTM_13N_WGS84_G2139_3D.wkt -r cubic -tr 1.0 1.0 merged_dsm.tif merged_dsm_reprojected_UTM_13N_WGS84_G2139.tif None """ @@ -88,16 +91,23 @@ def create_dsm(extent_geojson: str, # processing_extent.geojson src_wkt = f.read() POINTCLOUD_CRS = [src_wkt for _ in range(len(readers))] - + output_path = os.dirname(output_prefix) + prefix = os.path.basename(output_prefix) output_path = Path(output_path) output_path.mkdir(exist_ok=True) print(f"Number of readers: {len(readers)}") + dsm_fn_list = [] + dtm_fn_list = [] + intensity_fn_list = [] for i, reader in enumerate(readers): print(f"Processing reader #{i}") - dsm_file = output_path / f'dsm_tile_aoi_{str(i).zfill(4)}.tif' - dtm_file = output_path / f'dtm_tile_aoi_{str(i).zfill(4)}.tif' - intensity_file = output_path / f'intensity_tile_aoi_{str(i).zfill(4)}.tif' + dsm_file = output_path / f'{prefix}_dsm_tile_aoi_{str(i).zfill(4)}.tif' + dtm_file = output_path / f'{prefix}_dtm_tile_aoi_{str(i).zfill(4)}.tif' + intensity_file = output_path / f'{prefix}_intensity_tile_aoi_{str(i).zfill(4)}.tif' + dsm_fn_list.append(dsm_file) + dtm_fn_list.append(dtm_file) + intensity_fn_list.append(intensity_file) ## DSM creation block pipeline_dsm = {'pipeline':[reader]} @@ -122,7 +132,7 @@ def create_dsm(extent_geojson: str, # processing_extent.geojson gridmethod=GRID_METHOD, dimension='Z') pipeline_dsm['pipeline'] += pdal_pipeline_dsm pipeline_dsm['pipeline'] += dsm_stage - + # Save a copy of each pipeline dsm_pipeline_config_fn = os.path.join(output_path,f"pipeline_dsm_{str(i).zfill(4)}.json") with open(dsm_pipeline_config_fn, 'w') as f: @@ -200,3 +210,22 @@ def create_dsm(extent_geojson: str, # processing_extent.geojson f.write(json.dumps(pipeline_intensity)) pipeline_intensity = pdal.Pipeline(json.dumps(pipeline_intensity)) pipeline_intensity.execute() + + if mosaic: + print("*** Now creating raster composites ***") + dsm_mos_fn = f"{out_prefix}-DSM_mos.tif" + print(f"Creating DSM mosaic at {dsm_mos_fn}") + _ = dsm_functions.dem_mosaic(dsm_fn_list,dsm_mos_fn) + dtm_mos_fn = f"{out_prefix}-DTM_mos.tif" + print(f"Creating DTM mosaic at {dtm_mos_fn}") + _ = dsm_functions.dem_mosaic(dtm_fn_list,dtm_mos_fn) + intensity_mos_fn = f"{out_prefix}-Intensity_mos.tif" + print(f"Creating intensity raster mosaic at {intensity_mos_fn}") + _ = dsm_functions.dem_mosaic(intensity_fn_list,intensity_mos_fn) + if cleanup: + print("User selected to remove intermediate tile outputs") + for fn in dsm_fn_list + dtm_fn_list + intensity_fn_list: + os.remove(fn) + + + From 569bdf3771145781d82af02fbd0c618e3cf240db Mon Sep 17 00:00:00 2001 From: ShashankBice Date: Fri, 21 Feb 2025 18:49:54 +0000 Subject: [PATCH 02/11] rename and re-structre to allow proper install --- src/{grid_pc => lidar_tools}/__init__.py | 2 +- src/{grid_pc => lidar_tools}/__main__.py | 0 src/{grid_pc => lidar_tools}/cli.py | 0 src/{grid_pc => lidar_tools}/dsm_functions.py | 0 src/{grid_pc => lidar_tools}/filter_percentile.py | 0 src/{grid_pc => lidar_tools}/pdal_pipeline.py | 10 +++++----- 6 files changed, 6 insertions(+), 6 deletions(-) rename src/{grid_pc => lidar_tools}/__init__.py (76%) rename src/{grid_pc => lidar_tools}/__main__.py (100%) rename src/{grid_pc => lidar_tools}/cli.py (100%) rename src/{grid_pc => lidar_tools}/dsm_functions.py (100%) rename src/{grid_pc => lidar_tools}/filter_percentile.py (100%) rename src/{grid_pc => lidar_tools}/pdal_pipeline.py (97%) diff --git a/src/grid_pc/__init__.py b/src/lidar_tools/__init__.py similarity index 76% rename from src/grid_pc/__init__.py rename to src/lidar_tools/__init__.py index f5ed453..3c5b660 100644 --- a/src/grid_pc/__init__.py +++ b/src/lidar_tools/__init__.py @@ -1,6 +1,6 @@ from lidar_tools import filter_percentile from importlib.metadata import version as _version -__version__ = _version("grid_pc") +__version__ = _version("lidar_tools") __all__ = ["filter_percentile"] diff --git a/src/grid_pc/__main__.py b/src/lidar_tools/__main__.py similarity index 100% rename from src/grid_pc/__main__.py rename to src/lidar_tools/__main__.py diff --git a/src/grid_pc/cli.py b/src/lidar_tools/cli.py similarity index 100% rename from src/grid_pc/cli.py rename to src/lidar_tools/cli.py diff --git a/src/grid_pc/dsm_functions.py b/src/lidar_tools/dsm_functions.py similarity index 100% rename from src/grid_pc/dsm_functions.py rename to src/lidar_tools/dsm_functions.py diff --git a/src/grid_pc/filter_percentile.py b/src/lidar_tools/filter_percentile.py similarity index 100% rename from src/grid_pc/filter_percentile.py rename to src/lidar_tools/filter_percentile.py diff --git a/src/grid_pc/pdal_pipeline.py b/src/lidar_tools/pdal_pipeline.py similarity index 97% rename from src/grid_pc/pdal_pipeline.py rename to src/lidar_tools/pdal_pipeline.py index a94323a..dbf90ef 100644 --- a/src/grid_pc/pdal_pipeline.py +++ b/src/lidar_tools/pdal_pipeline.py @@ -2,7 +2,7 @@ Generate a DSM from input polygon """ import os -from grid_pc import dsm_functions +from lidar_tools import dsm_functions import pdal from shapely.geometry import Polygon import geopandas as gpd @@ -34,11 +34,11 @@ def create_dsm(extent_geojson: str, # processing_extent.geojson - source_wkt: str | None, # SRS_CRS.wkt target_wkt: str, # UTM_13N_WGS84_G2139_3D.wkt output_prefix: str, # CO_ALS_proc/CO_3DEP_ALS - mosaic: bool | True, - cleanup: bool | True) -> None: + source_wkt: str = None, # SRS_CRS.wkt + mosaic: bool = True, + cleanup: bool = True) -> None: """ Create a Digital Surface Model (DSM) from a given extent and point cloud data. This function divides a region of interest into tiles and generates a DSM geotiff from EPT point clouds for each tile using PDAL @@ -91,7 +91,7 @@ def create_dsm(extent_geojson: str, # processing_extent.geojson src_wkt = f.read() POINTCLOUD_CRS = [src_wkt for _ in range(len(readers))] - output_path = os.dirname(output_prefix) + output_path = os.path.dirname(output_prefix) prefix = os.path.basename(output_prefix) output_path = Path(output_path) output_path.mkdir(exist_ok=True) From 6a823e73ab529e9f1a2aef961b74ec173afcd43c Mon Sep 17 00:00:00 2001 From: ShashankBice Date: Fri, 21 Feb 2025 20:14:15 +0000 Subject: [PATCH 03/11] small aoi to check functions --- notebooks/very_small_aoi.geojson | 36 ++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) create mode 100644 notebooks/very_small_aoi.geojson diff --git a/notebooks/very_small_aoi.geojson b/notebooks/very_small_aoi.geojson new file mode 100644 index 0000000..ad253b6 --- /dev/null +++ b/notebooks/very_small_aoi.geojson @@ -0,0 +1,36 @@ +{ + "type": "FeatureCollection", + "features": [ + { + "type": "Feature", + "properties": {}, + "geometry": { + "coordinates": [ + [ + [ + -106.50967466214205, + 38.80653868377519 + ], + [ + -106.50967466214205, + 38.78329966412687 + ], + [ + -106.4710996229505, + 38.78329966412687 + ], + [ + -106.4710996229505, + 38.80653868377519 + ], + [ + -106.50967466214205, + 38.80653868377519 + ] + ] + ], + "type": "Polygon" + } + } + ] +} From 9be01652d02e07734950b8fc235b8014b82c257d Mon Sep 17 00:00:00 2001 From: ShashankBice Date: Fri, 21 Feb 2025 20:19:27 +0000 Subject: [PATCH 04/11] fix typo --- src/lidar_tools/pdal_pipeline.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/lidar_tools/pdal_pipeline.py b/src/lidar_tools/pdal_pipeline.py index dbf90ef..0807840 100644 --- a/src/lidar_tools/pdal_pipeline.py +++ b/src/lidar_tools/pdal_pipeline.py @@ -213,13 +213,13 @@ def create_dsm(extent_geojson: str, # processing_extent.geojson if mosaic: print("*** Now creating raster composites ***") - dsm_mos_fn = f"{out_prefix}-DSM_mos.tif" + dsm_mos_fn = f"{output_prefix}-DSM_mos.tif" print(f"Creating DSM mosaic at {dsm_mos_fn}") _ = dsm_functions.dem_mosaic(dsm_fn_list,dsm_mos_fn) - dtm_mos_fn = f"{out_prefix}-DTM_mos.tif" + dtm_mos_fn = f"{output_prefix}-DTM_mos.tif" print(f"Creating DTM mosaic at {dtm_mos_fn}") _ = dsm_functions.dem_mosaic(dtm_fn_list,dtm_mos_fn) - intensity_mos_fn = f"{out_prefix}-Intensity_mos.tif" + intensity_mos_fn = f"{output_prefix}-intensity_mos.tif" print(f"Creating intensity raster mosaic at {intensity_mos_fn}") _ = dsm_functions.dem_mosaic(intensity_fn_list,intensity_mos_fn) if cleanup: From d41de0aa466e7e050eb2c1a3fa8b0198315f10ed Mon Sep 17 00:00:00 2001 From: ShashankBice Date: Fri, 21 Feb 2025 20:30:24 +0000 Subject: [PATCH 05/11] syntax fix --- src/lidar_tools/dsm_functions.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/lidar_tools/dsm_functions.py b/src/lidar_tools/dsm_functions.py index b9a27fb..b882fd3 100644 --- a/src/lidar_tools/dsm_functions.py +++ b/src/lidar_tools/dsm_functions.py @@ -7,7 +7,6 @@ import geopandas as gpd import requests import subprocess -import ast from shutil import which def return_readers(input_aoi, @@ -230,7 +229,7 @@ def dem_mosaic(img_list,outfn,tr=None,tsrs=None,stats=None,tile_size=None,extent if stats: dem_mosaic_opt.extend(['--{}'.format(stats)]) - if (tr is not None) & (ast.literal_eval(tr) is not None): + if tr: dem_mosaic_opt.extend(['--tr', str(tr)]) if tsrs: dem_mosaic_opt.extend(['--t_srs', tsrs]) From 410376b8f775b2403bfbe78b1e0072508f51b29a Mon Sep 17 00:00:00 2001 From: ShashankBice Date: Fri, 21 Feb 2025 20:33:56 +0000 Subject: [PATCH 06/11] remove posix path for command line --- src/lidar_tools/pdal_pipeline.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/lidar_tools/pdal_pipeline.py b/src/lidar_tools/pdal_pipeline.py index 0807840..509c3fb 100644 --- a/src/lidar_tools/pdal_pipeline.py +++ b/src/lidar_tools/pdal_pipeline.py @@ -105,9 +105,9 @@ def create_dsm(extent_geojson: str, # processing_extent.geojson dsm_file = output_path / f'{prefix}_dsm_tile_aoi_{str(i).zfill(4)}.tif' dtm_file = output_path / f'{prefix}_dtm_tile_aoi_{str(i).zfill(4)}.tif' intensity_file = output_path / f'{prefix}_intensity_tile_aoi_{str(i).zfill(4)}.tif' - dsm_fn_list.append(dsm_file) - dtm_fn_list.append(dtm_file) - intensity_fn_list.append(intensity_file) + dsm_fn_list.append(dsm_file.as_posix()) + dtm_fn_list.append(dtm_file.as_posix()) + intensity_fn_list.append(intensity_file.as_posix()) ## DSM creation block pipeline_dsm = {'pipeline':[reader]} From 9d139b4f6b3364e12deb1ba67ad1eadb8cdb7f60 Mon Sep 17 00:00:00 2001 From: ShashankBice Date: Fri, 21 Feb 2025 21:12:09 +0000 Subject: [PATCH 07/11] update tile buffer value to 5 m --- src/lidar_tools/dsm_functions.py | 3 ++- src/lidar_tools/pdal_pipeline.py | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/lidar_tools/dsm_functions.py b/src/lidar_tools/dsm_functions.py index b882fd3..d4f4582 100644 --- a/src/lidar_tools/dsm_functions.py +++ b/src/lidar_tools/dsm_functions.py @@ -14,11 +14,12 @@ def return_readers(input_aoi, pointcloud_resolution = 1, n_rows = 5, n_cols=5, - buffer_value=0): + buffer_value=5): """ This method takes a raster file and finds overlapping 3DEP data. It then returns a series of readers corresponding to non overlapping areas that can be used as part of further PDAL processing pipelines The method also returns the CRS specified i + The default buffer value is 5, and in units of m """ xmin, ymin, xmax, ymax = input_aoi.bounds x_step = (xmax - xmin) / n_cols diff --git a/src/lidar_tools/pdal_pipeline.py b/src/lidar_tools/pdal_pipeline.py index 509c3fb..2a98f92 100644 --- a/src/lidar_tools/pdal_pipeline.py +++ b/src/lidar_tools/pdal_pipeline.py @@ -83,7 +83,7 @@ def create_dsm(extent_geojson: str, # processing_extent.geojson # Specfying a buffer_value > 0 will generate overlapping DEM tiles, resulting in a seamless # final mosaicked DEM readers, POINTCLOUD_CRS = dsm_functions.return_readers(input_aoi, input_crs, - pointcloud_resolution = 1, n_rows=5, n_cols=5, buffer_value=0) + pointcloud_resolution = 1, n_rows=5, n_cols=5, buffer_value=5) # NOTE: if source_wkt is passed, override POINTCLOUD_CRSs if source_wkt: From 15ecf74aa1b787c56bc0d116057a88c47a50e9d8 Mon Sep 17 00:00:00 2001 From: ShashankBice Date: Fri, 21 Feb 2025 21:36:58 +0000 Subject: [PATCH 08/11] update usage after restructre --- README.md | 31 ++++++++++++++++++++++++++++++- 1 file changed, 30 insertions(+), 1 deletion(-) diff --git a/README.md b/README.md index 09dd1c4..a8df468 100644 --- a/README.md +++ b/README.md @@ -17,7 +17,36 @@ We recommend using [pixi](https://pixi.sh/latest/) to install a locked software git clone https://github.com/uw-cryo/lidar_tools.git cd lidar_tools pixi shell -pdal_pipeline notebooks/processing_extent.geojson notebooks/SRS_CRS.wkt notebooks/UTM_13N_WGS84_G2139_3D.wkt /tmp/dem +pdal_pipeline notebooks/very_small_aoi.geojson notebooks/UTM_13N_WGS84_G2139_3D.wkt ../CO_processing/CO_3DEP_ALS --source-wkt notebooks/SRS_CRS.wkt +``` + +```console +> pdal_pipeline --help + + Usage: pdal_pipeline [OPTIONS] EXTENT_GEOJSON TARGET_WKT OUTPUT_PREFIX + + Create a Digital Surface Model (DSM) from a given extent and point cloud data. This function divides a region of interest into tiles and generates a DSM geotiff from EPT point clouds for each + tile using PDAL + Parameters ---------- extent_geojson : str Path to the GeoJSON file defining the processing extent. source_wkt : str or None Path to the WKT file defining the source coordinate reference + system (CRS). If None, the CRS from the point cloud file is used. target_wkt : str Path to the WKT file defining the target coordinate reference system (CRS). output_prefix : str prefix + with directory name and filename prefix for the project (e.g., CO_ALS_proc/CO_3DEP_ALS) mosaic : bool Mosaic the output tiles using a weighted average algorithm cleanup: bool If true, + remove the intermediate tif files for the output tiles Returns ------- None + Notes ----- After running this function, reproject the final DEM with the following commands: 1. gdalwarp -s_srs SRS_CRS.wkt -t_srs UTM_13N_WGS84_G2139_3D.wkt -r cubic -tr 1.0 1.0 merged_dsm.tif + merged_dsm_reprojected_UTM_13N_WGS84_G2139.tif None + +╭─ Arguments ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────╮ +│ * extent_geojson TEXT [default: None] [required] │ +│ * target_wkt TEXT [default: None] [required] │ +│ * output_prefix TEXT [default: None] [required] │ +╰──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────╯ +╭─ Options ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────╮ +│ --source-wkt TEXT [default: None] │ +│ --mosaic --no-mosaic [default: mosaic] │ +│ --cleanup --no-cleanup [default: cleanup] │ +│ --install-completion Install completion for the current shell. │ +│ --show-completion Show completion for the current shell, to copy it or customize the installation. │ +│ --help Show this message and exit. │ +╰──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────╯ ``` ## Installation From ad23fb25002410fba2730407bf3166331f897fa6 Mon Sep 17 00:00:00 2001 From: ShashankBice Date: Fri, 21 Feb 2025 21:49:03 +0000 Subject: [PATCH 09/11] fix buffer update step --- src/lidar_tools/dsm_functions.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/lidar_tools/dsm_functions.py b/src/lidar_tools/dsm_functions.py index d4f4582..de23d4c 100644 --- a/src/lidar_tools/dsm_functions.py +++ b/src/lidar_tools/dsm_functions.py @@ -41,7 +41,9 @@ def return_readers(input_aoi, aoi_3857 = Polygon.from_bounds(*src_bounds_transformed_3857) print(aoi.bounds, src_bounds_transformed_3857) if buffer_value: - aoi_3857.buffer(buffer_value) + aoi_3857 = aoi_3857.buffer(buffer_value) + print(f"The tile polygon is buffered by {buffer_value:.2f} m") + gdf = gpd.read_file('https://raw.githubusercontent.com/hobuinc/usgs-lidar/master/boundaries/resources.geojson').set_crs(4326) # in the eventuality that the above URL breaks, we store a local copy From 4841e2ae386e02a268581e98314651583f26c449 Mon Sep 17 00:00:00 2001 From: ShashankBice Date: Fri, 21 Feb 2025 21:57:50 +0000 Subject: [PATCH 10/11] update message for buffering --- src/lidar_tools/dsm_functions.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/lidar_tools/dsm_functions.py b/src/lidar_tools/dsm_functions.py index de23d4c..e882796 100644 --- a/src/lidar_tools/dsm_functions.py +++ b/src/lidar_tools/dsm_functions.py @@ -42,7 +42,7 @@ def return_readers(input_aoi, print(aoi.bounds, src_bounds_transformed_3857) if buffer_value: aoi_3857 = aoi_3857.buffer(buffer_value) - print(f"The tile polygon is buffered by {buffer_value:.2f} m") + print(f"The tile polygon will be buffered by {buffer_value:.2f} m") gdf = gpd.read_file('https://raw.githubusercontent.com/hobuinc/usgs-lidar/master/boundaries/resources.geojson').set_crs(4326) From 91fcaa872a78998b142d039bfa8a981bcb0a170f Mon Sep 17 00:00:00 2001 From: Shashank Bhushan Date: Mon, 24 Feb 2025 12:50:30 -0500 Subject: [PATCH 11/11] Update src/lidar_tools/pdal_pipeline.py Co-authored-by: Scott Henderson --- src/lidar_tools/pdal_pipeline.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/lidar_tools/pdal_pipeline.py b/src/lidar_tools/pdal_pipeline.py index 2a98f92..b98017e 100644 --- a/src/lidar_tools/pdal_pipeline.py +++ b/src/lidar_tools/pdal_pipeline.py @@ -102,9 +102,9 @@ def create_dsm(extent_geojson: str, # processing_extent.geojson intensity_fn_list = [] for i, reader in enumerate(readers): print(f"Processing reader #{i}") - dsm_file = output_path / f'{prefix}_dsm_tile_aoi_{str(i).zfill(4)}.tif' - dtm_file = output_path / f'{prefix}_dtm_tile_aoi_{str(i).zfill(4)}.tif' - intensity_file = output_path / f'{prefix}_intensity_tile_aoi_{str(i).zfill(4)}.tif' + dsm_file = f'{output_path}/{prefix}_dsm_tile_aoi_{str(i).zfill(4)}.tif' + dtm_file = f'{output_path}/{prefix}_dtm_tile_aoi_{str(i).zfill(4)}.tif' + intensity_file = f'{output_path}/{prefix}_intensity_tile_aoi_{str(i).zfill(4)}.tif' dsm_fn_list.append(dsm_file.as_posix()) dtm_fn_list.append(dtm_file.as_posix()) intensity_fn_list.append(intensity_file.as_posix())