uw-cryo · ShashankBice · Feb 24, 2025 · Feb 21, 2025 · Feb 21, 2025 · Feb 21, 2025
diff --git 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

diff --git a/notebooks/very_small_aoi.geojson 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"
+      }
+    }
+  ]
+}
diff --git a/src/grid_pc/__init__.py → src/lidar_tools/__init__.py b/src/grid_pc/__init__.py → 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 → src/lidar_tools/__main__.py b/src/grid_pc/__main__.py → src/lidar_tools/__main__.py
diff --git a/src/grid_pc/cli.py → src/lidar_tools/cli.py b/src/grid_pc/cli.py → src/lidar_tools/cli.py
diff --git a/src/grid_pc/dsm_functions.py → src/lidar_tools/dsm_functions.py b/src/grid_pc/dsm_functions.py → src/lidar_tools/dsm_functions.py
@@ -6,17 +6,20 @@
 from shapely.geometry import Polygon
 import geopandas as gpd
 import requests
+import subprocess
+from shutil import which
 
 def return_readers(input_aoi,
                    src_crs,
                    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
@@ -38,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 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)
             # in the eventuality that the above URL breaks, we store a local copy
@@ -198,3 +203,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:
+        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)
 dependencies = [ 
 dependencies = [ 
+    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/filter_percentile.py → src/lidar_tools/filter_percentile.py b/src/grid_pc/filter_percentile.py → src/lidar_tools/filter_percentile.py
diff --git a/src/grid_pc/pdal_pipeline.py → src/lidar_tools/pdal_pipeline.py b/src/grid_pc/pdal_pipeline.py → 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,10 +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_path: str, # /tmp/dem/
-               ) -> None:
+               output_prefix: str, # CO_ALS_proc/CO_3DEP_ALS
+               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
@@ -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
     """
 
@@ -80,24 +83,31 @@ 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:
         with open(source_wkt, 'r') as f:
             src_wkt = f.read()
         POINTCLOUD_CRS = [src_wkt for _ in range(len(readers))]
 
-
+    output_path = os.path.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 = 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())
 
         ## 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"{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"{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"{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:
+            print("User selected to remove intermediate tile outputs")
+            for fn in dsm_fn_list + dtm_fn_list + intensity_fn_list:
+                os.remove(fn)
+
+
+