notebooks/07_3DEP_Generate_Colorized_PointClouds.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2f30b87d",
   "metadata": {},
   "source": [
    "# Generate colorized (RGB) point clouds using USGS 3D Elevation Program (3DEP) lidar data and National Agriculture Imagery Program (NAIP) Imagery"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5e0f9db",
   "metadata": {},
   "source": [
    "## Table of Contents\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Purpose\" data-toc-modified-id=\"Purpose-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Purpose</a></span></li><li><span><a href=\"#Setup\" data-toc-modified-id=\"Setup-2\"><span class=\"toc-item-num\">2&nbsp;&nbsp;</span>Setup</a></span><ul class=\"toc-item\">\n",
    "<li><span><a href=\"#Installation-Options\" data-toc-modified-id=\"Installation-Options-2.1\"><span class=\"toc-item-num\">2.1&nbsp;&nbsp;</span>Installation Options</a></span><ul class=\"toc-item\"><li><span><a href=\"#Option-1:-Install-and-run-on-Google-Colaboratory\" data-toc-modified-id=\"Option-1:-Install-and-run-on-Google-Colaboratory-2.1.1\"><span class=\"toc-item-num\">2.1.1&nbsp;&nbsp;</span>Option 1: Install and run on Google Colaboratory</a></span></li><li><span><a href=\"#Option-2:-Install-and-run-on-local-file-system\" data-toc-modified-id=\"Option-2:-Install-and-run-on-local-file-system-2.1.2\"><span class=\"toc-item-num\">2.1.2&nbsp;&nbsp;</span>Option 2: Install and run on local file system</a></span></ul></li><li><span><a href=\"#Library-Imports\" data-toc-modified-id=\"Library-Imports-2.2\"><span class=\"toc-item-num\">2.2&nbsp;&nbsp;</span>Library Imports</a></span></li><li><span><a href=\"#Define-Functions\" data-toc-modified-id=\"Define-Functions-2.3\"><span class=\"toc-item-num\">2.3&nbsp;&nbsp;</span>Define Functions</a></span></ul></li><li><span><a href=\"#Data-Access-and-Processing\" data-toc-modified-id=\"Data-Access-and-Processing-3\"><span class=\"toc-item-num\">3&nbsp;&nbsp;</span>Data Access and Processing</a></span><ul class=\"toc-item\"><li><span><a href=\"#Get-3DEP-Dataset-Boundary-Polygons\" data-toc-modified-id=\"Get-3DEP-Dataset-Boundary-Polygons-3.1\"><span class=\"toc-item-num\">3.1&nbsp;&nbsp;</span>Get 3DEP Dataset Boundary Polygons</a></span></li><li><span><a href=\"#Create-Interactive-Ipyleaflet-Map-And-Define-AOI\" data-toc-modified-id=\"Create-Interactive-Ipyleaflet-Map-And-Define-AOI-3.2\"><span class=\"toc-item-num\">3.2&nbsp;&nbsp;</span>Create Interactive Ipyleaflet Map and Define AOI</a></span></li><li><span><a href=\"#Find-3DEP-Polygon(s)-Intersecting-AOI\" data-toc-modified-id=\"Find-3DEP-Polygon(s)-Intersecting-AOI-3.3\"><span class=\"toc-item-num\">3.3&nbsp;&nbsp;</span>Find 3DEP Polygon(s) Intersecting AOI</a></span></li><li><span><a href=\"#Specify-Point-Cloud-Resolution\" data-toc-modified-id=\"Specify-Point-Cloud-Resolution-3.4\"><span class=\"toc-item-num\">3.4&nbsp;&nbsp;</span>Specify Point Cloud Resolution</a></span></li><li><span><a href=\"#Query-National-Agriculture-Imagery-Program-Service-for-Imagery\" data-toc-modified-id=\"Query-National-Agriculture-Imagery-Program-Service-for-Imagery-3.5\"><span class=\"toc-item-num\">3.5&nbsp;&nbsp;</span>Query National Agriculture Imagery Program Service for Imagery</a></span></li><li><span><a href=\"#Construct-and-Exectute-PDAL-Pipeline\" data-toc-modified-id=\"|-3.6\"><span class=\"toc-item-num\">3.6&nbsp;&nbsp;</span>Construct and Exectute PDAL Pipeline</a></span></ul></li><li><span><a href=\"#Conclusion\" data-toc-modified-id=\"Conclusion-4\"><span class=\"toc-item-num\">4&nbsp;&nbsp;</span>Conclusion</a></span></li><li><span><a href=\"#Resources\" data-toc-modified-id=\"Resources-5\"><span class=\"toc-item-num\">5&nbsp;&nbsp;</span>Resources</a></span></li><li><span><a href=\"#Funding-Keywords-Citation\" data-toc-modified-id=\"Funding-Keywords-Citation-6\"><span class=\"toc-item-num\">6&nbsp;</span>Funding, Keywords, and Citation</a></span></li></ul></div>\n",
    "\n",
    "<a name=\"Authors\"></a>\n",
    "## Authors\n",
    "Cole Speed<sup>1</sup>, Matthew Beckley<sup>1</sup>, Christopher Crosby<sup>1</sup>, Viswanath Nandigam<sup>2</sup>\n",
    "\n",
    "<sup>1</sup>UNAVCO, Inc.; <sup>2</sup>San Diego Supercomputer Center\n",
    "\n",
    "<a name=\"Purpose\"></a>\n",
    "    \n",
    "## Purpose\n",
    "\n",
    "The United States Geological Survey's 3-D Elevation Program (3DEP) is collecting high-quality light detection and ranging (lidar) data over the conterminous U.S., Hawaii, and the U.S territories. Data acquisition is ongoing, with over 1800 3DEP datasets consisting of more than 42 trillion points covering an area of greater than 6.5 million $km^{2}$ already available for use. Data are hosted and publically available in Entwine Point Tile (EPT) format in an<a href=\"https://registry.opendata.aws/usgs-lidar/\"> Amazon Web Services (AWS) S3 public bucket</a>.\n",
    "\n",
    "3DEP point cloud data can be accessed through several user interfaces including the <a href=\"https://portal.opentopography.org/datasets\">OpenTopography Data Portal</a> and the <a href=\"https://prd-tnm.s3.amazonaws.com/LidarExplorer/index.html#/process\">USGS Lidar Explorer</a>. In many cases, however, users may prefer to access and produce derivative products from the 3DEP lidar data programmatically, particularly in the case of more advanced users seeking to access and process immense swaths of lidar point cloud data on local, commercial cloud, or high performance compute infrastructure. However, to date tools and/or workflows for programmatic access and processing of the 3DEP data remain underdeveloped.\n",
    "\n",
    "This Jupyter Notebook is part of a suite developed by OpenTopography in collaboration with the USGS National Geospatial Program, in an effort to leverage available APIs and cloud resources to enhance access and usability of 3DEP data and products for the geospatial community. Importantly, these notebooks are designed to be accessible by user's with limited experience in the Python programming language and as educational tools for those looking to learn about accessing, processing, and visualizing cloud-hosted data with Python. These notebooks use open-source Python libraries, including the Point Data Abstraction Libary (PDAL) and Geospatial Data Abstraction Library (GDAL) for all point cloud and raster processing. \n",
    "\n",
    "This and the other use-case specific Jupyter Notebooks developed as part of this effort for programmatic 3DEP access are available in a <a href=\"https://github.com/OpenTopography/OT_3DEP_Workflows\"> Github repository </a> and may be run locally or on the <a href=\"https://colab.research.google.com/\">Google Colaboratory</a> cloud platform.\n",
    "\n",
    "#### Specific features of this notebook\n",
    "\n",
    "1. Users may (1) select an area of interest (AOI) directly in an interactive map within this notebook or (2) import a shapefile defining their area of interest.\n",
    "\n",
    "2. API request to <a href=\"https://registry.opendata.aws/usgs-lidar/\"> Amazon Web Services (AWS) EPT (Entwine Point Tile) S3 bucket</a> returns 3DEP point cloud data within user-defined AOI for subsequent processing. \n",
    "\n",
    "3. Flexible and customizable PDAL pipelines are available for specifying point cloud resolution, filtering, reclassifying, and reprojecting points.\n",
    " \n",
    "4. API request to NAIP Imagery Service for imagery within the user-defined AOI.\n",
    "\n",
    "5. 3DEP lidar points are assigned RGB values from NAIP imagery\n",
    "\n",
    "5. Colorized point cloud data are saved to local directory or on Google Drive (if user uses Google Colaboratory option, see below.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04d9b4c8",
   "metadata": {},
   "source": [
    "<a name=\"Setup\"></a>\n",
    "## Setup\n",
    "\n",
    "<a name=\"Installation-Options\"></a>\n",
    "### Installation Options\n",
    "There are two options for performing the workflow steps outlined below. **Option 1** is our suggested method for simplicity, as building a virtual environment with the required dependencies on the user's local file system can be challenging based on the user's experience level with Python and <a href=\"https://www.anaconda.com/\"> Anaconda</a>.\n",
    "\n",
    "1. **Option 1**: Launch the interactive Jupyter notebook on Google Colaboratory.\n",
    "    - Does not require creation of a virtual environment or installation on local file system.\n",
    "    - Requires Google account and access to personal Google Drive folder.\n",
    "    - Data download limits will be dependent on user's available Google Drive storage. \n",
    "    - If you wish to run this notebook in Google Colaboratory click the 'Open in Colab' badge below.\n",
    "    - <font color='red'> **Note:** RAM allotment for Google Colaboratory may not be sufficient for very accessing large point clouds (>100,000,000 points) </font> \n",
    "    <br/><br/>\n",
    "2. **Option 2**: Download this Jupyter notebook (.ipynb file) to your local file system.\n",
    "    - Create a virtual conda environment containing the required dependencies.\n",
    "    - Run Juypter notebook on local machine.\n",
    "    - Data download limits and computation speed will be dependent on user's hardware.\n",
    "    - <font color='red'> **Note:** We suggest running notebooks locally on high-RAM hardware for jobs exceeding ~100,000,000 points </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d1db532",
   "metadata": {},
   "source": [
    "<a name=\"Option-1:-Install-and-run-on-Google-Colaboratory\"></a>\n",
    "### Option 1: Install and run on Google Colaboratory\n",
    "For ease-of-use, it is suggested to launch and execute these notebooks on <a href=\"https://colab.research.google.com/\">Google Colaboratory</a> (Colab, for short), Google's Cloud Platform. Dependencies will be installed on a virtual machine on Google's cloud servers and the code will be executed directly in your browser! A major benefit of this is that you will have direct access to Google's high-end CPU/GPUs and will not have to install any dependencies locally. All deliverables will be saved to your personal Google Drive. To experiment and run one of the below Jupyter Notebooks on Google Colab click the \"Open in Colab\" badge below.\n",
    "\n",
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/OpenTopography/OT_3DEP_Workflows/blob/main/notebooks/07_3DEP_Generate_Colorized_PointClouds.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "ca3dadb4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This cell only excecutes if you're running on Colab. Installation process takes 2-3 minutes.\n",
    "import os, sys\n",
    "if 'google.colab' in sys.modules:\n",
    "    \n",
    "  # Mount Google Drive. You will be prompted to grant file I/O access to Drive.\n",
    "  from google.colab import drive \n",
    "  drive.mount('/gdrive/') # Mount Google Drive! \n",
    "\n",
    "  # Clone OpenTopography 3DEP Workflow Git Repository\n",
    "  !git clone https://github.com/OpenTopography/OT_3DEP_Workflows\n",
    "\n",
    "  #  Install the core dependencies (other than PDAL/GDAL) from requirements.txt\n",
    "  !pip install -r OT_3DEP_Workflows/requirements.txt\n",
    "\n",
    "  # Install Conda (necessary to install PDAL/GDAL)\n",
    "  !pip install -q condacolab\n",
    "  import condacolab\n",
    "  condacolab.install()\n",
    "\n",
    "  #kernel will restart. Install PDAL and GDAL with Mamba.\n",
    "  !mamba install -q python-pdal gdal\n",
    "  \n",
    "  # Runtime will restart automatically. Do not rerun above cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "6957ba1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This cell only excecutes if you're running on Colab.\n",
    "import os, sys\n",
    "if 'google.colab' in sys.modules:\n",
    "    # Colab requires proj_lib environment variable to be set manually.\n",
    "    os.environ['PROJ_LIB'] = '/usr/local/share/proj/'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "455d380e",
   "metadata": {},
   "source": [
    "**If using Option 1 (Google Colab), proceed to Library Imports**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ff85a4d",
   "metadata": {},
   "source": [
    "<a name=\"Option-2:-Install-and-run-on-local-file-system\"></a>\n",
    "### Option 2: Install and run on local file system\n",
    "\n",
    "If you would like to run the Jupyter Notebook on your local machine:\n",
    "\n",
    "Make a new directory on your local file system where the 3DEP Jupyter Notebooks (and all 3DEP data, if desired) will be saved. In this example case, the directory will be called `3DEP`.\n",
    "  \n",
    "    $ mkdir 3DEP\n",
    "\n",
    "Change into the new directory and `git clone` the Github repository containing the Jupyter Notebooks and other relevant files to your local file system.\n",
    "\n",
    "    $ cd 3DEP\n",
    "    \n",
    "    $ git clone https://github.com/OpenTopography/OT_3DEP_Workflows\n",
    "\n",
    "Anaconda is recommended for Python package installation and management. Package versions in Anaconda are managed by the package management system *conda*. Anaconda installers for MacOS/Linux/Windows can be downloaded from https://docs.anaconda.com/anaconda/install/. Follow the instructions to install Anaconda.\n",
    "\n",
    "After installing Anaconda, create a virtual environment with the required dependencies, using the `environment.yml` file contained in the cloned Github repo). Note: Exectuting the following command will automatically create the conda environment with name `3dep` and all of the required dependencies installed. If you would prefer a different name, replace `3dep` with another name in the following command:\n",
    "\n",
    "\t$ conda env create -n 3dep --file environment.yml\n",
    "\n",
    "Next, activate the conda environment with all of the necessary dependencies installed. \n",
    "\t\n",
    "\t$ conda activate 3dep\n",
    "    \n",
    "Now, launch the chosen Jupyter Notebook. If unsure how to launch a Notebook, refer to this guide (https://jupyter-notebook-beginner-guide.readthedocs.io/en/latest/execute.html). \n",
    "\n",
    "**You may now proceed to Library Imports**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94d5bf63",
   "metadata": {},
   "source": [
    "<a name=\"Library-Imports\"></a>\n",
    "### Library Imports\n",
    "\n",
    "After successfully completing the steps outlined in either **Option 1** or **Option 2**, we can now import the modules for use throughout the rest of the notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ac6ab591",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Packages (Local Machine Only)\n",
    "import copy\n",
    "import geopandas as gpd\n",
    "import ipyleaflet\n",
    "import ipywidgets as widgets\n",
    "import json\n",
    "import math\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import os\n",
    "from osgeo import gdal\n",
    "import pdal\n",
    "import pdb\n",
    "import pyproj\n",
    "import requests\n",
    "from shapely.geometry import shape, Point, Polygon\n",
    "from shapely.ops import transform\n",
    "import subprocess\n",
    "import sys"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11f865f7",
   "metadata": {},
   "source": [
    "Now, set the necessary GDAL configuration options."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9afcb1e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "gdal.SetConfigOption('GDAL_HTTP_UNSAFESSL', 'YES')\n",
    "gdal.SetConfigOption('GDAL_ENABLE_DEPRECATED_DRIVER_DODS', 'YES')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af8b1e62",
   "metadata": {},
   "source": [
    "<a name=\"Define-Functions\"></a>\n",
    "### Define Functions\n",
    "\n",
    "Several functions are provided in the cell below. These functions are necessary for successful execution of remainder of the notebook. Broadly, these functions provide the utility for the user to draw and area of interest (AOI) on an interactive map and construct the PDAL pipeline for getting the point cloud data from the Amazon Web Services EPT bucket, performing processing steps, producing DEMs, and saving the results. A description of the parameters are provided as docstrings in the function definitions.\n",
    "\n",
    "\n",
    "**These functions can be modified as the user sees fit; however, they are designed to work with a simple execution of the below cell.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "3b544699",
   "metadata": {},
   "outputs": [],
   "source": [
    "def proj_to_3857(poly, orig_crs):\n",
    "    \"\"\"\n",
    "    Function for reprojecting a polygon from a shapefile of any CRS to Web Mercator (EPSG: 3857).\n",
    "    The original polygon must have a CRS assigned.\n",
    "    \n",
    "    :param poly: shapely polygon for user area of interest (AOI)\n",
    "    :param orig_crs: the original CRS for the shapefile. It is stripped out during import_shapefile_to_shapely() method\n",
    "    \"\"\"\n",
    "    wgs84 = pyproj.CRS(\"EPSG:4326\")\n",
    "    web_mercator = pyproj.CRS(\"EPSG:3857\")\n",
    "    project_gcs = pyproj.Transformer.from_crs(orig_crs, wgs84, always_xy=True).transform\n",
    "    project_wm = pyproj.Transformer.from_crs(orig_crs, web_mercator, always_xy=True).transform\n",
    "    user_poly_proj4326 = transform(project_gcs, poly)\n",
    "    user_poly_proj3857 = transform(project_wm, poly)\n",
    "    return(user_poly_proj4326, user_poly_proj3857)\n",
    "\n",
    "def gcs_to_proj(bbox):\n",
    "    \"\"\"\n",
    "    Function for reprojecting polygon shapely object from geographic coordinates (EPSG:4326)\n",
    "    to Web Mercator (EPSG: 3857)).\n",
    "\n",
    "    :param poly: shapely polygon for user area of interest (AOI)\n",
    "    \"\"\"\n",
    "    wgs84 = pyproj.CRS(\"EPSG:4326\")\n",
    "    web_mercator = pyproj.CRS(\"EPSG:3857\")\n",
    "    project = pyproj.Transformer.from_crs(wgs84, web_mercator, always_xy=True).transform\n",
    "    user_poly_proj3857 = transform(project, bbox)\n",
    "    return (user_poly_proj3857)\n",
    "\n",
    "def import_shapefile_to_shapely(path):\n",
    "    \"\"\"\n",
    "    Conversion of shapefile to shapely object.\n",
    "    \n",
    "    :param path: location of shapefile on user's local file system\n",
    "    \"\"\"\n",
    "    shapefile_path = path\n",
    "    gdf = gpd.read_file(shapefile_path)\n",
    "    orig_crs = gdf.crs                   # this is the original CRS of the imported shapefile\n",
    "    user_shp = gdf.loc[0, 'geometry']\n",
    "    user_shp_epsg4326, user_shp_epsg3857 = proj_to_3857(user_shp, orig_crs)\n",
    "    user_AOI = [[user_shp_epsg4326, user_shp_epsg3857]]\n",
    "    return user_AOI\n",
    "\n",
    "def handle_draw(target, action, geo_json):\n",
    "    \"\"\"\n",
    "    Functionality to draw area of interest (AOI) on interactive ipyleaflet map.\n",
    "    \n",
    "    :param extent_epsg3857: polygon for user-defined AOI\n",
    "    :param usgs_3dep_dataset_name: name of 3DEP dataset which AOI overlaps\n",
    "    :param resolution: The desired resolution of the pointcloud based on the following definition:\n",
    "    \"\"\"\n",
    "        \n",
    "    geom = dict(geo_json['geometry'])\n",
    "    user_poly = shape(geom)\n",
    "    user_poly_proj3857 = gcs_to_proj(user_poly)\n",
    "    print('AOI is valid and has boundaries of ', user_poly_proj3857.bounds, 'Please proceed to the next cell.')\n",
    "    user_AOI.append((user_poly, user_poly_proj3857))  #for various reasons, we need user AOI in GCS and EPSG 3857\n",
    "    \n",
    "\n",
    "def build_pdal_pipeline(extent_epsg3857, usgs_3dep_dataset_names, pc_resolution, filterNoise = True,\n",
    "                        colorize = True, savePointCloud = True, outCRS = 3857, rasterName = 'naip.tif',\n",
    "                        pc_outName = 'pointcloud_colorized', pc_outType = 'laz'):\n",
    "    \n",
    "    \"\"\"\n",
    "    Build pdal pipeline for requesting, processing, and saving point cloud data. Each processing step is a 'stage' \n",
    "    in the final pdal pipeline. Each stages is appended to the 'pointcloud_pipeline' object to produce the final pipeline.\n",
    "    \n",
    "    Parameters:\n",
    "    :param extent_epsg3857: Polygon for user-defined AOI in Web Mercator projection (EPS:3857)Polygon is generated \n",
    "                            either through the 'handle_draw' methor or by inputing their own shapefile.\n",
    "    :param usgs_3dep_dataset_names: List of name of the 3DEP dataset(s) that the data will be obtained. This parameter is set \n",
    "                                determined through intersecttino of the 3DEP and AOI polys.\n",
    "    :param pc_resolution: The desired resolution of the pointcloud based on the following definition:\n",
    "        \n",
    "                        Source: https://pdal.io/stages/readers.ept.html#readers-ept\n",
    "                            A point resolution limit to select, expressed as a grid cell edge length. \n",
    "                            Units correspond to resource coordinate system units. For example, \n",
    "                            for a coordinate system expressed in meters, a resolution value of 0.1 \n",
    "                            will select points up to a ground resolution of 100 points per square meter.\n",
    "                            The resulting resolution may not be exactly this value: the minimum possible \n",
    "                            resolution that is at least as precise as the requested resolution will be selected. \n",
    "                            Therefore the result may be a bit more precise than requested.\n",
    "                            \n",
    "    :param filterNoise: Option to remove points from USGS Class 7 (Low Noise) and Class 18 (High Noise).\n",
    "    :param colorize: Option to assign RGB values from NAIP imagery to the 3DEP lidar points.\n",
    "    :param savePointCloud: Option to save (or not) the point cloud data. If savePointCloud == False, \n",
    "           the pc_outName and pc_outType parameters are not used and can be any value.\n",
    "    :param outCRS: Output coordinate reference systemt (CRS), specified by ESPG code (e.g., 3857 - Web Mercator)\n",
    "    :param rasterName: Name of the NAIP tile that corresponds to the user AOI.\n",
    "    :param pc_outName: Desired name of file on user's local file system. If savePointcloud = False, \n",
    "                  pc_outName can be in value.\n",
    "    :param pc_outType:  Desired file extension. Input must be either 'las' or 'laz'. If savePointcloud = False, \n",
    "                  pc_outName can be in value. If a different file type is requested,the user will get error.\n",
    "    :raise Exception: If user passes in argument that is not 'las' or 'laz'.\n",
    "    \"\"\"\n",
    "    readers = []\n",
    "    for name in usgs_3dep_dataset_names:\n",
    "        url = \"https://s3-us-west-2.amazonaws.com/usgs-lidar-public/{}/ept.json\".format(name)\n",
    "        reader = {\n",
    "            \"type\": \"readers.ept\",\n",
    "            \"filename\": str(url),\n",
    "            \"polygon\": str(extent_epsg3857),\n",
    "            \"requests\": 3,\n",
    "            \"resolution\": pc_resolution\n",
    "        }\n",
    "        readers.append(reader)\n",
    "\n",
    "    pointcloud_pipeline = {\n",
    "        \"pipeline\":\n",
    "            readers\n",
    "    }\n",
    "\n",
    "    if filterNoise == True:\n",
    "        filter_stage = {\n",
    "            \"type\": \"filters.range\",\n",
    "            \"limits\": \"Classification![7:7], Classification![18:18]\"\n",
    "        }\n",
    "\n",
    "        pointcloud_pipeline['pipeline'].append(filter_stage)\n",
    "\n",
    "    if colorize == True:\n",
    "        colorize_stage = {\n",
    "                \"type\": \"filters.colorization\",\n",
    "                \"raster\": rasterName\n",
    "            }\n",
    "        color_filter_stage = {\n",
    "                \"type\": \"filters.range\",\n",
    "                \"limits\": \"Red[1:]\"\n",
    "            }\n",
    "        color_writer_stage = {\n",
    "                \"type\": \"writers.las\",\n",
    "                \"compression\": \"laszip\",\n",
    "                \"minor_version\": \"2\",\n",
    "                \"dataformat_id\": \"3\",\n",
    "                \"filename\": pc_outName\n",
    "            }\n",
    "\n",
    "        pointcloud_pipeline['pipeline'].append(colorize_stage)\n",
    "        pointcloud_pipeline['pipeline']. append(color_filter_stage)\n",
    "        pointcloud_pipeline['pipeline'].append(color_writer_stage)\n",
    "\n",
    "    if savePointCloud == True:\n",
    "\n",
    "        if pc_outType == 'las':\n",
    "            savePC_stage = {\n",
    "                \"type\": \"writers.las\",\n",
    "                \"filename\": str(pc_outName) + '.' + str(pc_outType)\n",
    "            }\n",
    "        elif pc_outType == 'laz':\n",
    "            savePC_stage = {\n",
    "                \"type\": \"writers.las\",\n",
    "                \"compression\": \"laszip\",\n",
    "                \"filename\": str(pc_outName) + '.' + str(pc_outType)\n",
    "            }\n",
    "        else:\n",
    "            raise Exception(\"pc_outType must be 'las' or 'laz'.\")\n",
    "\n",
    "        pointcloud_pipeline['pipeline'].append(savePC_stage)\n",
    "        \n",
    "        print(\"pdal pipeline built successfully...\")\n",
    "\n",
    "    return pointcloud_pipeline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0485d62e",
   "metadata": {},
   "source": [
    "<a name=\"Data-Access-and-Processing\"></a>\n",
    "## Data Access and Processing\n",
    "Now that we have the required modules imported and functions defined, we can proceed with defining our area of interest (AOI), accessing/processing the 3DEP data from the Amazon Web Services EPT bucket. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e850f0f4",
   "metadata": {},
   "source": [
    "<a name=\"Get-3DEP-Dataset-Boundary-Polygons\"></a>\n",
    "### Get 3DEP Dataset Boundary Polygons  \n",
    "First, we will get the 3DEP dataset polygon boundaries to see if 3DEP data is currently available for our area of interest. An up-to-date version of the currently available 3DEP dataset boundaries are maintained by Hobu Inc. in the usgs-lidar Github repository https://github.com/hobuinc/usgs-lidar/. \n",
    "\n",
    "In the following cell, we use an API request to get the boundaries from the repository and save a local copy of these boundaries in geojson format. Next, we create a geopandas dataframe object to easily access the names, url, and point count of each dataset, and we use the ```gcs_to_proj()``` function to project each 3DEP polygon geometry to Web Mercator projection (EPSG: 3857), which is the native projection of the 3DEP data in the AWS S3 bucket."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "47f8da16",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get GeoJSON file for 3DEP outlines from URL \n",
    "\n",
    "print(\"Requesting, loading, and projecting 3DEP dataset polygons...\")\n",
    "\n",
    "#request the boundaries from the Github repo and save locally.\n",
    "url = 'https://raw.githubusercontent.com/hobuinc/usgs-lidar/master/boundaries/resources.geojson'\n",
    "r = requests.get(url)\n",
    "with open('resources.geojson', 'w') as f:\n",
    "    f.write(r.content.decode(\"utf-8\"))\n",
    "\n",
    "with open('resources.geojson', 'r') as f:\n",
    "    geojsons_3DEP = json.load(f)\n",
    "    \n",
    "#make pandas dataframe and create pandas.Series objects for the names, urls, and number of points for each boundary.\n",
    "with open('resources.geojson', 'r') as f:\n",
    "    df = gpd.read_file(f)\n",
    "    names = df['name']\n",
    "    urls = df['url']\n",
    "    num_points = df['count']\n",
    "\n",
    "#project the boundaries to EPSG 3857 (necessary for API call to AWS for 3DEP data)\n",
    "projected_geoms = []\n",
    "for geometry in df['geometry']:\n",
    "        projected_geoms.append(gcs_to_proj(geometry))\n",
    "\n",
    "geometries_GCS = df['geometry']\n",
    "geometries_EPSG3857 = gpd.GeoSeries(projected_geoms)\n",
    "\n",
    "print('Done. 3DEP polygons downloaded and projected to Web Mercator (EPSG:3857)')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "349a5f46",
   "metadata": {},
   "source": [
    "<a name=\"Create-Interactive-Ipyleaflet-Map-And-Define-AOI\"></a>\n",
    "### Create Interactive Ipyleaflet Map and Define AOI\n",
    "Next, we will define our area of interest (AOI) using an ipyleaflet interactive map and the 3DEP dataset polygons loaded in the prior step.  There are currently two options for defining the AOI:\n",
    "\n",
    "1. **Option 1** - User loads their own shapefile (.shp) to define the AOI. If using Google Colab, shapefile should be uploaded to the Google Drive folder or within the runtime in the browser.\n",
    "\n",
    "- To use **option 1**, In the following cell,  the system path to the .shp should be included between the quotes (e.g., ```shapefile = '/path/to/shapefile.shp'```). The other associated files (.cpg, .dbf., .prj, .shx) must be present within the same folder.  If using Google Colab, the  path should be to the file in the user's Google Drive (i.e., /content/drive/shapefile.shp). If running locally, the path shoudl be to the location of the shapefile on the user's local file system.\n",
    "\n",
    "2. **Option 2** - User draws a polygon on an ipyleaflet interactive map to define the AOI.\n",
    "\n",
    " - To use **option 2**, the path should be left blank (```shapefile = ''```), as it is by default.\n",
    " \n",
    "**The cell must be run either way.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "4f801c40",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Enter shapefile path, if applicable. Example: shapefile_path = '/path/to/shapefile.shp'.\n",
    "# Otherwise leave as shapefile_path = ''\n",
    "# Run this cell either way, or next cell will not run appropriately.\n",
    "shapefile_path = ''"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cc48aa47",
   "metadata": {},
   "source": [
    "If **option 1** is chosen, running the next cell will load the shapefile and print `shapefile loaded. proceed to next cell`.\n",
    "\n",
    "If **option 2** is chosen, running the next cell will produce an interactive map, make an ipyleaflet layer out of the 3DEP dataset polygon geojson, and allow the user to draw an AOI. When the map is rendered, use the +/- symbols to zoom in and out. The currently available 3DEP datasets are shown in transparent green polygons. The current functionality of the notebook allows the user AOI to straddle multiple 3DEP datsets. In this case, points from both datasets will be queried and merged into the resulting las/laz file.\n",
    "\n",
    "**Important Note 1:** Use either the pentagon or rectangle shaped buttons on the left sidebar to define the AOI. The pentagon-shaped button will allow the user to draw a polygon of any shape, and the rectangle-shaped will result in a rectangular AOI. Either shape will work for the boundary parameter in the point cloud request, and either can be used to define the AOI.\n",
    "\n",
    "**Important Note 2:** If the user would like to draw a new AOI, the cell below must be re-executed each and every time. If multiple AOIs are drawn on the same interactive map without re-executing the cell, subsequent cells will not function appropriately. The cell must be run each and every time the user would like to define a new AOI.\n",
    "\n",
    "**Important Note 3:** There are currently no checks implemented to keep the user from defining an AOI encompassing billions of lidar points. However, if such an AOI is defined, the compute time and hard drive space required will be substantial. Depending on your personal computer specifications, massive point clouds may cause issues with available RAM. A later cell will provide an estimation of the total point count for the AOI, but it is advisable to keep the AOI ~5000 hectares or less. The area is shown on the interactive map when drawing the AOI.\n",
    "\n",
    "**Important Note 4 (Google Colaboratory Users Only):** Every Google Colab instance has a RAM allotment of 12GB. This should be plenty to perform the desired tasks, assuming that an AOI of reasonble size is defined (see Important Note 3 for what is meant by 'reasonable size'). Testing of maximum possible point cloud size allowable on Google Colab instances has been minimal(other than to determine that point clouds up to several hundred million points are possible). If you perform such tests, please let us know what you find!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "5a77d6cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This cell plots a map and allows the user to make an AOI\n",
    "\n",
    "m = ipyleaflet.Map(\n",
    "    basemap=ipyleaflet.basemaps.Esri.WorldTopoMap,\n",
    "    center=(37, -100),\n",
    "    zoom=3.5,\n",
    "    crs=ipyleaflet.projections.EPSG3857\n",
    "    )\n",
    "\n",
    "geo_json_3DEP = ipyleaflet.GeoJSON(data=geojsons_3DEP, style = {'color': 'green', 'opacity':1, \n",
    "                                       'weight':1, 'fillOpacity':0.1})\n",
    "\n",
    "dc = ipyleaflet.DrawControl(\n",
    "    polygon={\"shapeOptions\": {\"color\": \"blue\"}},\n",
    "    rectangle={\"shapeOptions\": {\"color\": \"blue\"}},\n",
    "    polyline={},\n",
    "    circle={},\n",
    "    marker={},\n",
    "    circlemarker={},\n",
    ")\n",
    "\n",
    "print('Select an Area of Interest using the tools on the left side of the map.')\n",
    "user_AOI = []\n",
    "dc.on_draw(handle_draw)\n",
    "m.add_control(dc)\n",
    "m.add_layer(geo_json_3DEP)  #add 3DEP polygons GeoJSON\n",
    "display(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a54e112a",
   "metadata": {},
   "source": [
    "<a name=\"Find-3DEP-Polygon(s)-Intersecting-AOI\"></a>\n",
    "### Find 3DEP Polygon(s) Intersecting AOI\n",
    "Now that the user-specified AOI is defined, the following cell will determine the intersecting 3DEP dataset names and show the corresponding polygons on an interactive map.  `intersecting_polys` will be a list of the intersecting 3DEP dataset name(s), boundary(ies) in EPSG: 4326, boundary(ies) in EPSG: 3857, url(s), and the number of points in the entire 3DEP dataset(s). The dataset names will be used in the API request to the AWS EPT S3 bucket. A ratio of the total number of points and the area of the user-defined AOI will be used to estimate the total points within the AOI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "ebc096f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "AOI_GCS = user_AOI[-1][0]\n",
    "AOI_EPSG3857 = user_AOI[-1][1]\n",
    "\n",
    "intersecting_polys = []\n",
    "\n",
    "for i,geom in enumerate(geometries_EPSG3857):\n",
    "    if geom.intersects(AOI_EPSG3857):\n",
    "        intersecting_polys.append((names[i], geometries_GCS[i], geometries_EPSG3857[i], urls[i], num_points[i]))\n",
    "        \n",
    "print(intersecting_polys)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7e7b485",
   "metadata": {},
   "source": [
    "<a name=\"Specify-Point-Cloud-Resolution\"></a>\n",
    "### Specify Point Cloud Resolution\n",
    "Executing the next cell will show the AOI, the relevant 3DEP dataset(s) on another interactive map, and the option to specify point cloud resolution. \n",
    "\n",
    "Importantly, after the map is rendered, the user must define the desired 'point cloud resolution' using the radio buttons below the map. An estimation of the total number of lidar points within the bounding box is provided based on the area of the AOI and the total number of lidar points in the 3DEP dataset(s). Selecting the \"Full\" option will return all points in the quad (this number can be quite large, depending on the size of the AOI). Selecting any of the other options for resolution will return the appropriate number of points to ensure at least one lidar point per Nth meter (where N is the chosen resolution). The Entwine Point Tile (EPT) file format utilizes an octree structure for the point cloud, and the chosen resolution defines how deep in the octree to request points to obtain the specified resolution. This depth, and total number points varies drastically based on a number of parameters including local topography and vegetation. Therefore, the 'resolution' paramater and the total point count do not scale linearly. In other words, specifying a resolution of 2m will likely return far less than half of the number of points returned with 'full' resolution. The estimate of the full poin total provided is not exact, but should give the user some idea of how many points to expect the resultant point cloud to contain. \n",
    "\n",
    "Select the appropriate radio button below the map to specify `pointcloud_resolution`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "057e9762",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find AOI center for plotting purposes\n",
    "centroid =  list(AOI_GCS.centroid.coords)[0]\n",
    "\n",
    "#make ipyleaflet map\n",
    "m = ipyleaflet.Map(\n",
    "    basemap=ipyleaflet.basemaps.Esri.WorldTopoMap,\n",
    "    center=(centroid[1],centroid[0]),\n",
    "    zoom=12,\n",
    "    )\n",
    "\n",
    "#add intersecting 3DEP polygon(s) to the map\n",
    "wlayer_3DEP_list = []\n",
    "usgs_3dep_datasets = []\n",
    "number_pts_est = []\n",
    "for i, poly in enumerate(intersecting_polys):\n",
    "    wlayer_3DEP = ipyleaflet.WKTLayer(\n",
    "        wkt_string=poly[1].wkt, \n",
    "        style={\"color\": \"green\"})\n",
    "    \n",
    "    m.add_layer(wlayer_3DEP)\n",
    "    wlayer_3DEP_list.append(wlayer_3DEP)\n",
    "    usgs_3dep_datasets.append(poly[0])\n",
    "    \n",
    "    #estimate total points using ratio of area and point count\n",
    "    number_pts_est.append((int((AOI_EPSG3857.area/poly[2].area)*(poly[4]))))\n",
    "\n",
    "\n",
    "#make ipyleaflet layers from the AOI and add to map\n",
    "wlayer_user = ipyleaflet.WKTLayer(\n",
    "    wkt_string=AOI_GCS.boundary.wkt,\n",
    "    style={\"color\": \"blue\"}\n",
    ")\n",
    "\n",
    "AOI_EPSG3857_wkt = AOI_EPSG3857.wkt\n",
    "m.add_layer(wlayer_user)\n",
    "\n",
    "\n",
    "#sum the estimates of the number of points from each 3DEP dataset within the AOI\n",
    "num_pts_est = sum(number_pts_est)\n",
    "\n",
    "#Plot map and specify desired point cloud resolution using a widget\n",
    "user_resolution = widgets.RadioButtons(\n",
    "    options=[\n",
    "        (f'Full - All ~{int(math.ceil(num_pts_est/1e6)*1e6):,} points', 1.0),\n",
    "        (f'High - 2m resolution', 2.0),\n",
    "        (f'Mid  - 5m resolution', 5.0),\n",
    "        (f'Low  - 10m resolution', 10.0)\n",
    "    ],\n",
    "    layout={'width': 'max-content'},\n",
    "    disabled = False,\n",
    ")\n",
    "\n",
    "display(m)\n",
    "print(f'Your AOI at full resolution will include approximately {int(math.ceil(num_pts_est/1e6)*1e6):,} points. Select desired point cloud resolution.')\n",
    "widgets.VBox(\n",
    "    [user_resolution]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5dad0558",
   "metadata": {},
   "source": [
    "<font color='red'>**Note**: Lidar point clouds can get *very* large, *very* fast, and the point counts provided above are an estimate. Selecting the `Full` or `High` option when using Google Colaboratory (12GB RAM allocation) may cause the runtime to fail when attempting to access data exceeding ~50-100 million points. If full resolution is desired, we recommend running this notebook locally on hardware with more RAM. Keep this in mind when deciding the AOI size and point cloud resolution above!</font>\n",
    "\n",
    "The AOI bounding box, the relevant 3DEP dataset name(s), and the desired point cloud resolution are now defined. We can proceed with the API request to the AWS EPT bucket, processing, visualizing, and saving the data."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a051d065",
   "metadata": {},
   "source": [
    "<a name=\"Query-National-Agriculture-Imagery-Program-Service-for-Imagery\"></a>\n",
    "## Query National Agriculture Imagery Program Service for Imagery\n",
    "The National Agriculture Imagery Program (NAIP) acquires aerial imagery at a resolution of 1-meter ground sample distance (GSD) for the United States during the agricultural growing season. The imagery is available as a map service hosted in a number of locations including the <a href=\"https://services.nationalmap.gov/arcgis/rest/services/USGSNAIPPlus/MapServer\">USGS NAIP Plus service</a>. The data are originally tiled in 3.75' quadrangles.  We can query this service using the user-defined AOI to get a subset of the NAIP imagery corresponding to our point cloud location. \n",
    "\n",
    "First we must build our call to the NAIP service REST API."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "c96dfcaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define spatial extent\n",
    "#'coords' must be formatted as a list. I.e., coords = [minlon, minlat, maxlon, maxlat]\n",
    "minlon = AOI_GCS.bounds[0]\n",
    "minlat = AOI_GCS.bounds[1]\n",
    "maxlon = AOI_GCS.bounds[2]\n",
    "maxlat = AOI_GCS.bounds[3]\n",
    "\n",
    "bb = (minlon, minlat, maxlon, maxlat)\n",
    "bbox = ','.join(map(str, bb))         # put into CSV string for API access\n",
    "\n",
    "#Now define the USGS NAIP Plus service\n",
    "url = 'https://tnmaccess.nationalmap.gov/api/v1/products'\n",
    "product = 'USDA National Agriculture Imagery Program (NAIP)'\n",
    "extent = '3.75 x 3.75 minute'\n",
    "fmt = 'JPEG2000'\n",
    "version = 1\n",
    "\n",
    "#make a 'params' variable for the REST API call\n",
    "params = dict(datasets=product, bbox=bbox, prodFormats=fmt, outputFormat='JSON', version=version)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c4fe6c8",
   "metadata": {},
   "source": [
    "Now exectute the call to the USGS NAIP Plus REST API. There are several possibilities here if the following step does not execute properly. There is a print statement implemented that should provide some indication of which it is. \n",
    "\n",
    "1. If you receive the error `Error with Service Call` or `Error loading JSON output`, this most likely indicates that the region you selected does not have NAIP imagery associated with it. \n",
    "\n",
    "2. If you recieve the error `{'message': 'Endpoint request timed out'}`, this likely indicates the NAIP service is down. To check, click this link (https://stats.uptimerobot.com/gxzRZFARLZ/783928857). If you are greeted with the message `ScienceBase is down.`, the service is temporarily down and it is not possible to query the service at this time. \n",
    "\n",
    "If the cell executes successfully, information about the intersecting NAIP tile(s) will be printed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "d488a173",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Requesting NAIP imagery from the USGS NAIP Plus service...\")\n",
    "\n",
    "# Try the request with the params we have defined.\n",
    "try:\n",
    "    r = requests.get(url, params=params)\n",
    "except:\n",
    "    print('Error with Service Call')    #this would indicate that imagery is not available for the AOI\n",
    "    pdb.set_trace()\n",
    "\n",
    "# load API JSON output into a variable\n",
    "try:\n",
    "    data = json.loads(r.content)\n",
    "    print(data)\n",
    "except:\n",
    "    print('Error loading JSON output')\n",
    "    pdb.set_trace()\n",
    "\n",
    "# extract just the records from the JSON\n",
    "items = data['items']\n",
    "\n",
    "# extract the download URLS..\n",
    "URLS = []\n",
    "for f in items:\n",
    "    download_url = '/vsicurl/'+f['downloadURL']\n",
    "    URLS.append(download_url)\n",
    "\n",
    "# write URLs to a file (the AOI will often straddle multiple 3.75' tiles, this gets the names of all of them)\n",
    "# ---------------------------------\n",
    "outf = 'AOI_tiff_URLs.txt'\n",
    "out_obj = open(outf, 'w')\n",
    "\n",
    "for file in URLS:\n",
    "    out_obj.write(file + '\\n')\n",
    "out_obj.close()\n",
    "\n",
    "# ---------------------------------\n",
    "\n",
    "# Call GDALbuildvrt via system call.  Build a temporary VRT from the list\n",
    "# of files.  This vrt can be discarded after job is completed.\n",
    "vrt = \"tmp.vrt\"\n",
    "cmd1 = ['time gdalbuildvrt ' + vrt + ' -input_file_list ' + outf]\n",
    "\n",
    "# if running on colab this will execute the shell command\n",
    "if 'google.colab' in sys.modules:\n",
    "    !{cmd1[0]}\n",
    "# if not running on colab this will execute the shell command\n",
    "else:\n",
    "    p1 = subprocess.run(cmd1, shell=True)\n",
    "    if (p1.returncode == 1) or (p1.stderr is not None):\n",
    "        print('Problem building the VRT')\n",
    "        pdb.set_trace()\n",
    "\n",
    "# Call gdal_translate to subset the imagery (via the VRT). \n",
    "#Modify the outf name to something specific (e.g., boulderCO_naip.tif)\n",
    "outf = 'naip.tif'\n",
    "\n",
    "# projwin form: <ulx> <uly> <lrx> <lry>\n",
    "projwin_box = str(minlon) + ' ' + str(maxlat) + ' ' + str(maxlon) + ' ' + str(minlat)\n",
    "cmd2 = ['time gdal_translate ' + vrt + ' -projwin ' + projwin_box + ' -projwin_srs EPSG:4326 -of GTiff -co \"TILED=YES\" ' + outf]\n",
    "\n",
    "# if running on colab this will execute the shell command\n",
    "if 'google.colab' in sys.modules:\n",
    "    !{cmd2[0]}\n",
    "# if not running on colab this will execute the shell command\n",
    "else:\n",
    "    p2 = subprocess.run(cmd2, shell=True)\n",
    "    if (p2.returncode == 1) or (p2.stderr is not None):\n",
    "        print('Problem running gdal_translate')\n",
    "        pdb.set_trace()\n",
    "print(\"NAIP for your AOI has been saved\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d53a746c",
   "metadata": {},
   "source": [
    "If the above cell ran successfully, you now have the NAIP imagery for the AOI extent and can proceed with requesting the 3DEP data and colorizing the lidar points by the pixel values."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ddde0e6a",
   "metadata": {},
   "source": [
    "<a name=\"Construct-and-Exectute-PDAL-Pipeline\"></a>\n",
    "### Construct and Exectute PDAL Pipeline for Point Cloud Data\n",
    "The Point Data Abstraction Library (PDAL) is an open-source tool for translating and manipulating point cloud data. PDAL pipelines are useful ways of processing and manipulating point cloud data and creating derivative products. Pipelines comprise one or more stages that are read and executed in order on the point cloud dataset(s). \n",
    "\n",
    "The PDAL pipeline is constructed using the ```build_pdal_pipeline()``` function, and will construct the appropriate pipeline for the user's specifications (defined as function arguments). Executing this pipeline will make the API request, perform processing on the point cloud data (chosen by user) and provide the final result on the user's file system of Google Drive (Google Colab).  \n",
    "\n",
    "Paramaters (for more detailed descriptions of parameters, see <a href=\"#Define-Functions\" data-toc-modified-id=\"Define-Functions-6.3\">the function definitions</a>, above: <br>\n",
    "```AOI_EPSG3857_wkt```: the user-defined area of interest (AOI)<br>\n",
    "```usgs_3dep_dataset_names```: the intersecting 3DEP dataset names<br>\n",
    "```pointcloud_resolution```: point cloud resolution (1m, 2m, 5m, 10m)<br>\n",
    "```filterNoise```: remove the points of Class 7 (low noise) and Class 18 (high noise);<br>\n",
    "```colorize```: add rgb values from NAIP imagery to the points within each pixel.<br>\n",
    "```savePointCloud```: specify if point cloud data should be saved to local file system<br>\n",
    "```outCRS```: specify the coordinate reference system (CRS, in EPSG) of the output dataset.<br>\n",
    "```rasterName```: name of the NAIP raster file on local file system.\n",
    "```pc_outName```: name of point cloud on local file system<br>\n",
    "```pc_outType```: file type, las or laz (laszip compression). Options are 'las' or 'laz'<br>\n",
    "\n",
    "**Important Note 1: The ```AOI_EPSG3857_wkt```, ```usgs_3dep_datasets```, and ```pointcloud_resolution``` arguments are already defined. These should not be modified.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "dce80caf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Do not modify AOI_EPSG3857_wkt, usgs_3dep_datasets, or pointcloud_resolution\n",
    "# Modify the optional arguments to fit user need.\n",
    "# Change outCRS to EPSG code of desired coordinate reference system (Default is EPSG:3857 - Web Mercator Projection)\n",
    "# Change pc_outname to descriptive name and pc_outType to 'las' or 'laz'.\n",
    "\n",
    "pointcloud_resolution = user_resolution.value \n",
    "pc_pipeline = build_pdal_pipeline(AOI_EPSG3857_wkt, usgs_3dep_datasets, pointcloud_resolution, filterNoise = True,\n",
    "                        colorize = True, savePointCloud = True, outCRS = 3857, rasterName = 'naip.tif',\n",
    "                        pc_outName = 'pointcloud_colorized', pc_outType = 'laz')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c71aa286",
   "metadata": {},
   "source": [
    "The PDAL pipeline is now constructed. Running the the PDAL Python bindings function ```pdal.Pipeline()``` creates the pdal.Pipeline object from a json-ized version of the pointcloud pipeline we created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "1eb3c6cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "pc_pipeline = pdal.Pipeline(json.dumps(pc_pipeline))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62f44ed9",
   "metadata": {},
   "source": [
    "The cell below will execute the pc_pipeline object, which will make the API request, performing processing, and save the point cloud (if `savePointCloud == True`) at the specified location, name, and extension.\n",
    "\n",
    "The pipeline is executed in streaming mode, which significantly speeds up the process and cuts down on the required RAM. The `%%time` magic command will return the total computation time. The final output is the total number of points returned. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ee6acc37",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "pc_pipeline.execute_streaming(chunk_size=1000000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f95c4ca5",
   "metadata": {},
   "source": [
    "If the above cell ran successfully, the colorized point cloud should be saved on the user's local file system. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "159e51c0",
   "metadata": {},
   "source": [
    "<a name=\"Conclusion\"></a>\n",
    "## Conclusion\n",
    "This Jupyter Notebook is designed to make accessing, processing, visualizing, and creating derivative products with 3DEP more straightforward. Specifically, this notebook provides a workflow for programmatically accessing the 3DEP lidar point cloud data using an API request to the cloud-hosted AWS S3 bucket, performing basic processing steps on the point cloud data, and colorizing the lidar points by the corresponding NAIP imagery RGB values.\n",
    "\n",
    "Feedback regarding the workflow, specific use cases, and any computational testing with the notebook are welcomed and encouraged!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4808794d",
   "metadata": {},
   "source": [
    "<a name=\"Resources\"></a>\n",
    "## Resources\n",
    "\n",
    "All OpenTopography USGS 3DEP scientific workflows in this collection:<br>\n",
    "\n",
    "1. [Generate and visualize DEMs (DTM and DSM) from USGS 3D Elevation Program (3DEP) lidar data for user-defined area of interest](https://github.com/OpenTopography/OT_3DEP_Workflows/blob/main/notebooks/01_3DEP_Generate_DEM_User_AOI.ipynb)[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/OpenTopography/OT_3DEP_Workflows/blob/main/notebooks/01_3DEP_Generate_DEM_User_AOI.ipynb) <br>\n",
    "\n",
    "2. [Generate and visualize DEMs (DTM and DSM) from USGS 3D Elevation Program (3DEP) lidar data for USGS 7.5’ Quadrangles](https://github.com/OpenTopography/OT_3DEP_Workflows/blob/main/notebooks/02_3DEP_Generate_DEM_USGS_7.5_Quadrangles.ipynb)[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/OpenTopography/OT_3DEP_Workflows/blob/main/notebooks/02_3DEP_Generate_DEM_USGS_7.5_Quadrangles.ipynb) <br>\n",
    "\n",
    "3. [Generate and visualize DEMs (DTM and DSM) from USGS 3D Elevation Program (3DEP) lidar data for USGS Hydrologic Units](https://github.com/OpenTopography/OT_3DEP_Workflows/blob/main/notebooks/03_3DEP_Generate_DEM_USGS_HUCs.ipynb)[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/OpenTopography/OT_3DEP_Workflows/blob/main/notebooks/03_3DEP_Generate_DEM_USGS_HUCs.ipynb) <br>\n",
    "\n",
    "4. [Generate and visualize DEMs (DTM and DSM) from USGS 3D Elevation Program (3DEP) lidar data for user-defined corridors](https://github.com/OpenTopography/OT_3DEP_Workflows/blob/main/notebooks/04_3DEP_Generate_DEM_Corridors.ipynb)[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/OpenTopography/OT_3DEP_Workflows/blob/main/notebooks/04_3DEP_Generate_DEM_Corridors.ipynb) <br>\n",
    "\n",
    "5. [Generate Canopy Height Model (CHM) using USGS 3D Elevation Program (3DEP) lidar data for user-defined area of interest](https://github.com/OpenTopography/OT_3DEP_Workflows/blob/main/notebooks/05_3DEP_Generate_Canopy_Height_Models_User_AOI.ipynb)[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/OpenTopography/OT_3DEP_Workflows/blob/main/notebooks/05_3DEP_Generate_Canopy_Height_Models_User_AOI.ipynb) <br>\n",
    "\n",
    "6. [Topographic Differencing using USGS 3D Elevation Program (3DEP) lidar data for user-defined area of interest](https://github.com/OpenTopography/OT_3DEP_Workflows/blob/main/notebooks/06_3DEP_Topographic_Differencing.ipynb)[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/OpenTopography/OT_3DEP_Workflows/blob/main/notebooks/06_3DEP_Topographic_Differencing.ipynb) <br>\n",
    "\n",
    "7. [Generate colorized (RGB) point clouds using USGS 3D Elevation Program (3DEP) lidar data and National Agriculture Imagery Program (NAIP) Imagery](https://github.com/OpenTopography/OT_3DEP_Workflows/blob/main/notebooks/07_3DEP_Generate_Colorized_PointClouds.ipynb)[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/OpenTopography/OT_3DEP_Workflows/blob/main/notebooks/07_3DEP_Generate_Colorized_PointClouds.ipynb)\n",
    "\n",
    "### Additional Resources\n",
    "\n",
    "- Access USGS 3DEP via the <a href=\"https://portal.opentopography.org/datasets\">OpenTopography</a> portal (Currently restricted to academics).\n",
    "\n",
    "- The USGS 3DEP Lidar Point Cloud Data are accessible in Entwine Point Tile (EPT) format from this <a href=\"https://registry.opendata.aws/usgs-lidar/\">Amazon Web Services S3 Bucket</a>. <br>\n",
    "\n",
    "- The USGS hydrologic unit boundaries are accessed via the <a href=\"https://hydro.nationalmap.gov/arcgis/rest/services/wbd/MapServer\">USGS Watershed Dataset Map Service</a>. <br>\n",
    "\n",
    "- The USGS 7.5' quadrangle boundaries are accessed via the <a href=\"https://carto.nationalmap.gov/arcgis/rest/services/map_indices/MapServer\"> USGS Map Indicies Service</a>. <br>\n",
    "\n",
    "- Documentation for open-source Python libararies used by these workflows include <a href=\"https://pdal.dev/en/latest/\">PDAL</a> and <a href=\"https://gdal.org/\">GDAL</a>."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19b1a6e8",
   "metadata": {},
   "source": [
    "<a name=\"Funding-Keywords-Citation\"></a>\n",
    "## Funding\n",
    "\n",
    "Funding for the creation and distribution of these Jupyter Notebook-based workflows was provided as by the <a href=\"https://www.usgs.gov/centers/community-for-data-integration-cdi\">USGS Community for Data Integration (CDI)</a> through the funded grant <a href=\"https://www.usgs.gov/centers/community-for-data-integration-cdi/science/enhancing-usability-3dep-data-and-web-services\"> *Enhancing usability of 3DEP data and web services with Jupyter notebooks*</a>. OpenTopography is supported by the National Science Foundation (NSF) under Award Numbers <a href=\"https://nsf.gov/awardsearch/showAward?AWD_ID=1948997\">1948997</a>, <a href=\"https://nsf.gov/awardsearch/showAward?AWD_ID=1948994\">1948994</a> & <a href =\"https://nsf.gov/awardsearch/showAward?AWD_ID=1948857\">1948857</a>.\n",
    "\n",
    "## Keywords\n",
    "\n",
    "keywords=[\"OpenTopography\",\"USGS\", \"CDI\", \"3DEP\", \"PDAL\", \"Colorize point cloud\"]\n",
    "\n",
    "## Citation\n",
    "\n",
    "To cite this notebook:  Speed, C., Beckley, M., Crosby, C., & Nandigam, V. (2022). Generate colorized (RGB) point clouds using USGS 3D Elevation Program (3DEP) lidar data and National Agriculture Imagery Program (NAIP) Imagery (Version v1.0). DOI: Accessed: MM/DD/YYYY"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}