notebooks/0a_Download_LIDAR_data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given a raster file, we want to query and download available USGS 3DEP data. The data should conform to the bounds, the grid size, and projection of the raster file. Options are provided to filter out noise, and generate DTM/DSM, and saving the associated point cloud. We also provide an option to filter out LIDAR returns that exceed a user specified percentile threshold."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PDAL imports\n",
    "import pdal\n",
    "\n",
    "# rasterio imports\n",
    "import rasterio\n",
    "from rasterio.warp import transform_bounds\n",
    "from rasterio import warp\n",
    "from rasterio.merge import merge\n",
    "\n",
    "# GIS imports\n",
    "from pyproj import CRS\n",
    "from shapely.geometry import Polygon\n",
    "\n",
    "# Dataframe imports\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "\n",
    "# Misc imports\n",
    "import json\n",
    "import requests\n",
    "from pathlib import Path\n",
    "import os\n",
    "\n",
    "# set path to local ASP install for DEM mosaicking\n",
    "os.environ['PATH'] = os.environ['PATH'] + os.pathsep + '/mnt/working/karthikv/StereoPipeline-3.5.0-alpha-2024-10-15-x86_64-Linux/bin'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def return_readers(filename, n_rows = 5, n_cols=5, buffer_value=0):\n",
    "    \"\"\"\n",
    "    This method takes a raster file and finds overlapping 3DEP data. It then returns a series of readers\n",
    "    corresponding to non overlapping areas that can be used as part of further PDAL processing pipelines\n",
    "    The method also returns the CRS specified i\n",
    "    \"\"\"\n",
    "    with rasterio.open(filename) as ds:\n",
    "        src_bounds = ds.bounds\n",
    "        src_crs = ds.crs\n",
    "        src_transform = ds.transform\n",
    "\n",
    "    # this assumes the values to be in meters\n",
    "    x_resolution, y_resolution = abs(src_transform[0]), src_transform[4]\n",
    "    \n",
    "    # the point cloud resolution will be determined by the \n",
    "    # coarsest resolution available in our raster data\n",
    "    pointcloud_resolution = max([x_resolution, y_resolution])\n",
    "\n",
    "    xmin, ymin, xmax, ymax = src_bounds\n",
    "    x_step = (xmax - xmin) / n_cols\n",
    "    y_step = (ymax - ymin) / n_rows\n",
    "\n",
    "    dst_crs = CRS.from_epsg(4326)\n",
    "\n",
    "    readers = []\n",
    "\n",
    "    for i in range(int(n_cols)):\n",
    "        for j in range(int(n_rows)):\n",
    "            aoi = Polygon.from_bounds(xmin+i*x_step, ymin+j*y_step, xmin+(i+1)*x_step, ymin+(j+1)*y_step)\n",
    "\n",
    "            src_bounds_transformed = transform_bounds(src_crs, dst_crs, *aoi.bounds)\n",
    "            aoi_4326 = Polygon.from_bounds(*src_bounds_transformed)\n",
    "\n",
    "            src_bounds_transformed_3857 = transform_bounds(src_crs, CRS.from_epsg(3857), *aoi.bounds)\n",
    "            aoi_3857 = Polygon.from_bounds(*src_bounds_transformed_3857)\n",
    "            if buffer_value:\n",
    "                aoi_3857.buffer(buffer_value)\n",
    "\n",
    "            gdf = gpd.read_file('https://raw.githubusercontent.com/hobuinc/usgs-lidar/master/boundaries/resources.geojson').set_crs(4326) \n",
    "            # in the eventuality that the above URL breaks, we store a local copy\n",
    "            # gdf = gpd.read_file('../data/shapefiles/resources.geojson').set_crs(4326)\n",
    "\n",
    "            for _, row in gdf.iterrows():\n",
    "                if row.geometry.intersects(aoi_4326):\n",
    "                    usgs_dataset_name = row['name']\n",
    "                    break\n",
    "\n",
    "            url = f\"https://s3-us-west-2.amazonaws.com/usgs-lidar-public/{usgs_dataset_name}/ept.json\"\n",
    "            reader = {\n",
    "            \"type\": \"readers.ept\",\n",
    "            \"filename\": url,\n",
    "            \"resolution\": pointcloud_resolution,\n",
    "            \"polygon\": str(aoi_3857.wkt),\n",
    "            }\n",
    "            \n",
    "            # SRS associated with the 3DEP dataset\n",
    "            response = requests.get(url)\n",
    "            data = response.json()\n",
    "            srs_wkt = data['srs']['wkt']\n",
    "            pointcloud_input_crs = CRS.from_wkt(srs_wkt)\n",
    "            \n",
    "            readers.append(reader)\n",
    "\n",
    "    return readers, pointcloud_input_crs\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function that returns a PDAL pipeline to create a pointcloud based on user flags\n",
    "def create_pdal_pipeline(filter_low_noise=True, filter_high_noise=True, \n",
    "                         filter_road=True, reset_classes=False, reclassify_ground=False,\n",
    "                         return_only_ground=False, percentile_filter=True, percentile_threshold=0.95,\n",
    "                         group_filter=None,\n",
    "                         reproject=True, save_pointcloud=False, \n",
    "                         pointcloud_file = 'pointcloud', input_crs=None,\n",
    "                         output_crs=None, output_type='laz'):\n",
    "    \n",
    "    assert abs(percentile_threshold) <= 1, \"Percentile threshold must be in range [0, 1]\"\n",
    "    assert output_type in ['las', 'laz'], \"Output type must be either 'las' or 'laz'\"\n",
    "    assert output_crs is not None, \"Argument 'output_crs' must be explicitly specified!\"\n",
    "\n",
    "    stage_filter_low_noise = {\n",
    "        \"type\":\"filters.range\",\n",
    "        \"limits\":\"Classification![7:7]\"\n",
    "    }\n",
    "    stage_filter_high_noise = {\n",
    "        \"type\":\"filters.range\",\n",
    "        \"limits\":\"Classification![18:18]\"\n",
    "    }\n",
    "    stage_filter_road = {\n",
    "        \"type\":\"filters.range\",\n",
    "        \"limits\":\"Classification![11:11]\"\n",
    "    }\n",
    "    stage_reset_classes = {\n",
    "        \"type\":\"filters.assign\",\n",
    "        \"value\":\"Classification = 0\"\n",
    "    }\n",
    "    stage_reclassify_ground = {\n",
    "        \"type\":\"filters.smrf\",\n",
    "        # added from pdal smrf documentation, in turn from Pingel, 2013\n",
    "        \"scalar\":1.2,\n",
    "        \"slope\":0.2,\n",
    "        \"threshold\":0.45,\n",
    "        \"window\":8.0\n",
    "    }\n",
    "    stage_group_filter = {\n",
    "        \"type\":\"filters.returns\",\n",
    "        \"groups\":group_filter\n",
    "    }\n",
    "    stage_percentile_filter =  {\n",
    "        \"type\":\"filters.python\",\n",
    "        \"script\":\"filter_percentile.py\",\n",
    "        \"pdalargs\": {\"percentile_threshold\":percentile_threshold},\n",
    "        \"function\":\"filter_percentile\",\n",
    "        \"module\":\"anything\"\n",
    "    }\n",
    "    stage_return_ground = {\n",
    "        \"type\":\"filters.range\",\n",
    "        \"limits\":\"Classification[2:2]\"\n",
    "    }\n",
    "\n",
    "    # stage_reprojection = {\n",
    "    #     \"type\":\"filters.reprojection\",\n",
    "    #     \"out_srs\":str(output_crs)\n",
    "    # }\n",
    "    # if input_crs is not None:\n",
    "    #     stage_reprojection[\"in_srs\"] = str(input_crs)\n",
    "    \n",
    "    stage_save_pointcloud_las = {\n",
    "        \"type\": \"writers.las\",\n",
    "        \"filename\": f\"{pointcloud_file}.las\"\n",
    "    }\n",
    "    stage_save_pointcloud_laz = {\n",
    "        \"type\": \"writers.las\",\n",
    "        \"compression\": \"true\",\n",
    "        \"minor_version\": \"2\",\n",
    "        \"dataformat_id\": \"0\",\n",
    "        \"filename\": f\"{pointcloud_file}.laz\"\n",
    "    }\n",
    "\n",
    "    # Build pipeline\n",
    "    pipeline = []\n",
    "\n",
    "    # resetting the original classifications resets \n",
    "    # all point classifications to 0 (Unclassified)    \n",
    "    if reset_classes:\n",
    "        pipeline.append(stage_reset_classes)\n",
    "        if reclassify_ground:\n",
    "            pipeline.append(stage_reclassify_ground)\n",
    "    else:\n",
    "        # we apply the percentile filter first as it \n",
    "        # classifies detected outliers as 'high noise'\n",
    "        if group_filter is not None:\n",
    "            pipeline.append(stage_group_filter)\n",
    "        if percentile_filter:\n",
    "            pipeline.append(stage_percentile_filter)\n",
    "        if filter_low_noise:\n",
    "            pipeline.append(stage_filter_low_noise)\n",
    "        if percentile_filter or filter_high_noise:\n",
    "            pipeline.append(stage_filter_high_noise)\n",
    "        if filter_road:\n",
    "            pipeline.append(stage_filter_road)\n",
    "\n",
    "\n",
    "    # For creating DTMs, we want to process only ground returns\n",
    "    if return_only_ground:\n",
    "        pipeline.append(stage_return_ground)\n",
    "    \n",
    "    # if reproject:\n",
    "    #     pipeline.append(stage_reprojection)\n",
    "    \n",
    "    # the pipeline can save the pointclouds to a separate file if needed\n",
    "    if save_pointcloud:\n",
    "        if output_type == 'laz':\n",
    "            pipeline.append(stage_save_pointcloud_laz)\n",
    "        else:\n",
    "            pipeline.append(stage_save_pointcloud_las)\n",
    "\n",
    "    return pipeline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function that returns a PDAL pipeline to create a DEM based on user flags\n",
    "def create_dem_stage(dem_filename='dem_output.tif', pointcloud_resolution=1.,\n",
    "                        gridmethod='idw', dimension='Z'):\n",
    "    dem_stage = {\n",
    "            \"type\":\"writers.gdal\",\n",
    "            \"filename\":dem_filename,\n",
    "            \"gdaldriver\":'GTiff',\n",
    "            \"nodata\":-9999,\n",
    "            \"output_type\":gridmethod,\n",
    "            \"resolution\":float(pointcloud_resolution),\n",
    "            \"gdalopts\":\"COMPRESS=LZW,TILED=YES,blockxsize=256,blockysize=256,COPY_SRC_OVERVIEWS=YES\"\n",
    "    }\n",
    "\n",
    "    if dimension == 'Z':\n",
    "        dem_stage.update({\n",
    "            'dimension': 'Z',\n",
    "            'where': 'Z>0'\n",
    "        })\n",
    "    else:\n",
    "        dem_stage.update({\n",
    "            'dimension':dimension\n",
    "        })\n",
    "    \n",
    "    return [dem_stage]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# raster file for which pointcloud is generated\n",
    "input_file = '/mnt/working/karthikv/DeepDEM/data/09232024/WV01_20150911_1020010042D39D00_1020010043455300/original_rasters/1020010042D39D00_left_ortho_1.0m_final_aligned_to_REFDEM.tif'\n",
    "\n",
    "# we use a user specified output srs\n",
    "crs_file = '/mnt/working/karthikv/DeepDEM/data/09232024/WV01_20150911_1020010042D39D00_1020010043455300/original_rasters/UTM_10N_WGS84_G2139_3D.wkt'\n",
    "with open(crs_file, 'r') as f:\n",
    "    OUTPUT_CRS = CRS.from_string(f.read())\n",
    "\n",
    "# The method returns pointcloud readers, as well as the pointcloud file CRS as a WKT string\n",
    "# Specfying a buffer_value > 0 will generate overlapping DEM tiles, resulting in a seamless\n",
    "# final mosaicked DEM\n",
    "readers, POINTCLOUD_CRS = return_readers(input_file, n_rows=4, n_cols=4, buffer_value=100)\n",
    "\n",
    "# Set pointcloud processing parameters \n",
    "FILTER_LOW_NOISE = True\n",
    "FILTER_HIGH_NOISE = True\n",
    "FILTER_ROAD = True\n",
    "RETURN_ONLY_GROUND = False # Set true for DTM\n",
    "RESET_CLASSES = False\n",
    "RECLASSIFY_GROUND = False\n",
    "PERCENTILE_FILTER = True # Set to True to apply percentile based filtering of Z values\n",
    "PERCENTILE_THRESHOLD = 0.95 # Percentile value to filter out noisy Z returns\n",
    "# REPROJECT = True\n",
    "SAVE_POINTCLOUD=False\n",
    "POINTCLOUD_RESOLUTION = 1\n",
    "OUTPUT_TYPE='laz'\n",
    "GRID_METHOD='idw'\n",
    "DIMENSION='Z' # can be set to options accepted by writers.gdal. Set to 'intensity' to return intensity rasters\n",
    "\n",
    "# The group filter can be used to isolate first, intermediate, and last returns\n",
    "# Setting this to \"first,only\" will utilize top of canopy returns\n",
    "# set to \"None\" to use all returns\n",
    "GROUP_FILTER=\"first,only\"\n",
    "\n",
    "output_path = Path('/mnt/working/karthikv/DeepDEM/data/common_files/tmp')\n",
    "output_path.mkdir(exist_ok=True)\n",
    "\n",
    "for i, reader in enumerate(readers):\n",
    "    dem_file = output_path / f'dem_tile_aoi_{str(i).zfill(4)}.tif'\n",
    "    pipeline = {'pipeline':[reader]}\n",
    "\n",
    "    pdal_pipeline = create_pdal_pipeline(\n",
    "        filter_low_noise=FILTER_LOW_NOISE,\n",
    "        filter_high_noise=FILTER_HIGH_NOISE,\n",
    "        filter_road=FILTER_ROAD,\n",
    "        reset_classes=RESET_CLASSES, reclassify_ground=RECLASSIFY_GROUND,\n",
    "        return_only_ground=RETURN_ONLY_GROUND, \n",
    "        percentile_filter=PERCENTILE_FILTER, percentile_threshold=PERCENTILE_THRESHOLD,\n",
    "        group_filter=GROUP_FILTER,\n",
    "        # reproject=REPROJECT,\n",
    "        save_pointcloud=SAVE_POINTCLOUD, \n",
    "        pointcloud_file='pointcloud', input_crs = POINTCLOUD_CRS,\n",
    "        output_crs=OUTPUT_CRS, output_type=OUTPUT_TYPE\n",
    "    )\n",
    "\n",
    "    dem_stage = create_dem_stage(dem_filename=str(dem_file), \n",
    "                                    pointcloud_resolution=POINTCLOUD_RESOLUTION, \n",
    "                                    gridmethod=GRID_METHOD, dimension=DIMENSION)\n",
    "\n",
    "    # apply interpolation to fill gaps when generating DTM\n",
    "    if RETURN_ONLY_GROUND:\n",
    "        dem_stage[0]['window_size'] = 4\n",
    "    \n",
    "    pipeline['pipeline'] += pdal_pipeline\n",
    "    pipeline['pipeline'] += dem_stage\n",
    "    pipeline = pdal.Pipeline(json.dumps(pipeline))\n",
    "    pipeline.execute()\n",
    "\n",
    "# merge rasters \n",
    "if RETURN_ONLY_GROUND:\n",
    "    merged_filename = output_path.parent / 'merged_dtm.tif'\n",
    "else:\n",
    "    merged_filename = output_path.parent / 'merged_dsm.tif'\n",
    "\n",
    "os.system(f\"dem_mosaic -o {str(merged_filename)} {str(output_path)}/*.tif\")\n",
    "\n",
    "# delete temporary files\n",
    "for file in list(output_path.glob('*.tif')):\n",
    "    file.unlink()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "deep_dem",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}