datasets/3dep-lidar/3dep-lidar-cog-example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4ce5f387-9fb9-4adb-8f4f-d4bf7f92db17",
   "metadata": {},
   "source": [
    "## Accessing 3DEP Lidar COG data with the Planetary Computer STAC API\n",
    "\n",
    "The Planetary Computer includes a [group of datasets](https://planetarycomputer.microsoft.com/dataset/group/3dep-lidar) derived from the USGS 3DEP Lidar program. The raw data are available as a [collection of COPC assets](https://planetarycomputer.microsoft.com/dataset/3dep-lidar-copc). In addition, various derived products like [Intensity](https://planetarycomputer.microsoft.com/dataset/3dep-lidar-intensity) and [Height Above Ground](https://planetarycomputer.microsoft.com/dataset/3dep-lidar-hag) are available.\n",
    "\n",
    "This notebook demonstrates working with the derived COG products. For more on working with the COPC data, see its [example notebook](https://planetarycomputer.microsoft.com/dataset/3dep-lidar-copc#Example-Notebook).\n",
    "\n",
    "### Environment setup\n",
    "\n",
    "This notebook works with or without an API key, but you will be given more permissive access to the data with an API key.\n",
    "The [Planetary Computer Hub](https://planetarycomputer.microsoft.com/compute) is pre-configured to use your API key."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1236bb08-b2b2-48f7-9556-a46f41b1ba46",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pystac_client\n",
    "import planetary_computer\n",
    "import rasterio\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "# Set the environment variable PC_SDK_SUBSCRIPTION_KEY, or set it here.\n",
    "# The Hub sets PC_SDK_SUBSCRIPTION_KEY automatically.\n",
    "# pc.settings.set_subscription_key(<YOUR API Key>)\n",
    "\n",
    "client = pystac_client.Client.open(\n",
    "    \"https://planetarycomputer.microsoft.com/api/stac/v1/\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21e91ca0-da85-428e-a6b0-dfa481988e54",
   "metadata": {},
   "source": [
    "We'll use the STAC API to search for items from these various collections over Washington D.C."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "336a24ec-1c7e-4f16-b27d-471b8439f55c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'3dep-lidar-returns': <Item id=USGS_LPC_VA_Fairfax_County_2018-returns-5m-2-1>,\n",
       " '3dep-lidar-pointsourceid': <Item id=USGS_LPC_VA_Fairfax_County_2018-pointsourceid-5m-2-1>,\n",
       " '3dep-lidar-intensity': <Item id=USGS_LPC_VA_Fairfax_County_2018-intensity-2m-5-3>,\n",
       " '3dep-lidar-hag': <Item id=USGS_LPC_VA_Fairfax_County_2018-hag-2m-5-3>,\n",
       " '3dep-lidar-dtm-native': <Item id=USGS_LPC_VA_Fairfax_County_2018-dtm_native-2m-5-3>,\n",
       " '3dep-lidar-dtm': <Item id=USGS_LPC_VA_Fairfax_County_2018-dtm-2m-5-3>,\n",
       " '3dep-lidar-dsm': <Item id=USGS_LPC_VA_Fairfax_County_2018-dsm-2m-5-3>,\n",
       " '3dep-lidar-classification': <Item id=USGS_LPC_VA_Fairfax_County_2018-classification-5m-2-1>}"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "collections = [\n",
    "    \"3dep-lidar-hag\",\n",
    "    \"3dep-lidar-dsm\",\n",
    "    \"3dep-lidar-pointsourceid\",\n",
    "    \"3dep-lidar-intensity\",\n",
    "    \"3dep-lidar-dtm\",\n",
    "    \"3dep-lidar-dtm-native\",\n",
    "    \"3dep-lidar-returns\",\n",
    "    \"3dep-lidar-classification\",\n",
    "]\n",
    "\n",
    "search = client.search(\n",
    "    collections=collections,\n",
    "    intersects={\n",
    "        \"type\": \"Point\",\n",
    "        \"coordinates\": [-77.10058811018344, 38.838335717896314],\n",
    "    },\n",
    "    datetime=\"2018\",\n",
    ")\n",
    "items = planetary_computer.sign(search.get_all_items())\n",
    "items = {x.collection_id: x for x in items}\n",
    "\n",
    "items"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6bf3d93-cfb0-412f-8d35-0476117fafa8",
   "metadata": {},
   "source": [
    "### Height Above Ground\n",
    "\n",
    "This COG type is generated using the Z dimension of the [COPC data](https://planetarycomputer.microsoft.com/dataset/3dep-lidar-copc) data and removes noise, water, and using [`pdal.filters.smrf`](https://pdal.io/stages/filters.smrf.html#filters-smrf) followed by [pdal.filters.hag_nn](https://pdal.io/stages/filters.hag_nn.html#filters-hag-nn)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9cf656cd-40d6-4513-a694-6c662976af5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = rasterio.open(items[\"3dep-lidar-hag\"].assets[\"data\"].href).read().squeeze()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9fcdd3a0-4453-4f10-99b1-f68c6e465503",
   "metadata": {},
   "source": [
    "This COG has a special colormap we can use to visualize its values correctly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7c6c835e-d744-429c-ac4e-677d42bb3999",
   "metadata": {},
   "outputs": [],
   "source": [
    "pairs = [\n",
    "    ((-900, 1), (0, 0, 0, 0)),\n",
    "    ((1, 2), (205, 224, 241, 255)),\n",
    "    ((2, 3), (175, 209, 231, 255)),\n",
    "    ((3, 4), (137, 190, 220, 255)),\n",
    "    ((4, 5), (96, 166, 210, 255)),\n",
    "    ((5, 6), (34, 114, 181, 255)),\n",
    "    ((6, 7), (10, 84, 158, 255)),\n",
    "    ((7, 100), (8, 48, 107, 255)),\n",
    "]\n",
    "\n",
    "\n",
    "intervals, colors = zip(*pairs)\n",
    "nodes = np.array([x[1] for x in intervals]).astype(float)\n",
    "nodes -= np.abs(nodes.min())\n",
    "nodes /= nodes.max()\n",
    "\n",
    "\n",
    "colors = [np.asarray(c) / 255 for c in colors]\n",
    "\n",
    "cmap = matplotlib.colors.LinearSegmentedColormap.from_list(\n",
    "    \"hag\", list(zip(nodes, colors))\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "29c880bf-b0e8-4be7-9622-98c4726e98b3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"https://ai4edatasetspublicassets.blob.core.windows.net/assets/notebook-output/datasets-3dep-lidar-3dep-lidar-cog-example.ipynb/5.png\"/>"
      ],
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(8, 8))\n",
    "\n",
    "ax.imshow(cmap(ds), cmap=cmap)\n",
    "ax.set_axis_off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b074d8dc-dcc5-41c5-9481-1ac06a92c7d3",
   "metadata": {},
   "source": [
    "### Intensity\n",
    "\n",
    "This collection is derived from the [USGS 3DEP COPC collection](https://planetarycomputer.microsoft.com/dataset/3dep-lidar-copc). It is a collection of Cloud Optimized GeoTIFFs representing the pulse return magnitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b6cd8021-4257-46e4-98fe-9f286445fe46",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"https://ai4edatasetspublicassets.blob.core.windows.net/assets/notebook-output/datasets-3dep-lidar-3dep-lidar-cog-example.ipynb/6.png\"/>"
      ],
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ds = rasterio.open(items[\"3dep-lidar-intensity\"].assets[\"data\"].href).read().squeeze()\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 8))\n",
    "\n",
    "ax.imshow(ds, cmap=\"gray\")\n",
    "ax.set_axis_off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a67920ba-1afc-400d-b937-5cd6b3269b25",
   "metadata": {},
   "source": [
    "### Returns\n",
    "\n",
    "This collection is derived from the [USGS 3DEP COPC collection](https://planetarycomputer.microsoft.com/dataset/3dep-lidar-copc). It is a collection of Cloud Optimized GeoTIFFs representing the number of returns for a given pulse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c6825644-d141-4204-baa3-477366c45dae",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"https://ai4edatasetspublicassets.blob.core.windows.net/assets/notebook-output/datasets-3dep-lidar-3dep-lidar-cog-example.ipynb/7.png\"/>"
      ],
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ds = rasterio.open(items[\"3dep-lidar-returns\"].assets[\"data\"].href).read().squeeze()\n",
    "np.putmask(ds, ds < 1, 0)\n",
    "np.putmask(ds, ds >= 7, 7)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 8))\n",
    "\n",
    "ax.imshow(ds)\n",
    "ax.set_axis_off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e085816e-afee-4d21-9daa-de2c5009db07",
   "metadata": {},
   "source": [
    "### Classification\n",
    "\n",
    "This collection is derived from the [USGS 3DEP COPC collection](https://planetarycomputer.microsoft.com/dataset/3dep-lidar-copc). It uses the [ASPRS](https://www.asprs.org/) (American Society for Photogrammetry and Remote Sensing) [Lidar point classification](https://desktop.arcgis.com/en/arcmap/latest/manage-data/las-dataset/lidar-point-classification.htm). See [LAS specification](https://www.ogc.org/standards/LAS) for details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a154eafe-dc47-472a-98f2-78ff0ee7563c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"https://ai4edatasetspublicassets.blob.core.windows.net/assets/notebook-output/datasets-3dep-lidar-3dep-lidar-cog-example.ipynb/8.png\"/>"
      ],
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ds = (\n",
    "    rasterio.open(items[\"3dep-lidar-classification\"].assets[\"data\"].href)\n",
    "    .read()\n",
    "    .squeeze()\n",
    ")\n",
    "ds = np.where(ds > 0, ds, np.nan)\n",
    "fig, ax = plt.subplots(figsize=(8, 8))\n",
    "\n",
    "ax.imshow(ds, cmap=\"tab20\")\n",
    "ax.set_axis_off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd132505-5e18-4116-8f2f-72ebded42e96",
   "metadata": {},
   "source": [
    "### DSM\n",
    "\n",
    "This collection is derived from the [USGS 3DEP COPC collection](https://planetarycomputer.microsoft.com/dataset/3dep-lidar-copc). It creates a Digital Surface Model (DSM) using [`pdal.filters.range`](https://pdal.io/stages/filters.range.html#filters-range) to output a collection of Cloud Optimized GeoTIFFs, removing all points that have been classified as noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4f7484c0-8ec7-4428-ae76-6fb91c13180c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"https://ai4edatasetspublicassets.blob.core.windows.net/assets/notebook-output/datasets-3dep-lidar-3dep-lidar-cog-example.ipynb/9.png\"/>"
      ],
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ds = rasterio.open(items[\"3dep-lidar-dsm\"].assets[\"data\"].href).read().squeeze()\n",
    "ds = np.where(ds > 0, ds, np.nan)\n",
    "fig, ax = plt.subplots(figsize=(8, 8))\n",
    "\n",
    "ax.imshow(ds, cmap=\"gray\")\n",
    "ax.set_axis_off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af9b40ef-da80-4971-a96e-41a7be05049b",
   "metadata": {},
   "source": [
    "### DTM\n",
    "\n",
    "This collection is derived from the [USGS 3DEP COPC collection](https://planetarycomputer.microsoft.com/dataset/3dep-lidar-copc). It creates a Digital Terrain Model (DTM) using [`pdal.filters.smrf`](https://pdal.io/stages/filters.smrf.html#filters-smrf) to output a collection of Cloud Optimized GeoTIFFs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "3afbef14-f6e7-4ab5-a7aa-20b2801c1914",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"https://ai4edatasetspublicassets.blob.core.windows.net/assets/notebook-output/datasets-3dep-lidar-3dep-lidar-cog-example.ipynb/10.png\"/>"
      ],
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ds = rasterio.open(items[\"3dep-lidar-dtm\"].assets[\"data\"].href).read().squeeze()\n",
    "ds = np.where(ds > 0, ds, np.nan)\n",
    "fig, ax = plt.subplots(figsize=(8, 8))\n",
    "\n",
    "ax.imshow(ds, cmap=\"gray\")\n",
    "ax.set_axis_off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2c781381-53a1-49e1-bf70-9f4047dc621d",
   "metadata": {},
   "source": [
    "### DTM Native\n",
    "\n",
    "This collection is derived from the [USGS 3DEP COPC collection](https://planetarycomputer.microsoft.com/dataset/3dep-lidar-copc). It creates a Digital Terrain Model (DTM) using the vendor provided (native) ground classification and [`pdal.filters.range`](https://pdal.io/stages/filters.range.html#filters-range) to output a collection of Cloud Optimized GeoTIFFs, removing all points that have been classified as noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "48919fd8-39d0-4edd-ae65-61ce829771b7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"https://ai4edatasetspublicassets.blob.core.windows.net/assets/notebook-output/datasets-3dep-lidar-3dep-lidar-cog-example.ipynb/11.png\"/>"
      ],
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ds = rasterio.open(items[\"3dep-lidar-dtm-native\"].assets[\"data\"].href).read().squeeze()\n",
    "ds = np.where(ds > 0, ds, np.nan)\n",
    "fig, ax = plt.subplots(figsize=(8, 8))\n",
    "\n",
    "ax.imshow(ds, cmap=\"gray\")\n",
    "ax.set_axis_off()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e48c2fe-e8c7-4214-b898-f140eecdb020",
   "metadata": {},
   "source": [
    "### Point Source ID\n",
    "\n",
    "This collection is derived from the [USGS 3DEP COPC collection](https://planetarycomputer.microsoft.com/dataset/3dep-lidar-copc). It is a collection of Cloud Optimized GeoTIFFs representing the file source ID from which the point originated. Zero indicates that the point originated in the current file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "ae507018-d725-482c-b90e-1cebd5bc7741",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"https://ai4edatasetspublicassets.blob.core.windows.net/assets/notebook-output/datasets-3dep-lidar-3dep-lidar-cog-example.ipynb/12.png\"/>"
      ],
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ds = (\n",
    "    rasterio.open(items[\"3dep-lidar-pointsourceid\"].assets[\"data\"].href)\n",
    "    .read()\n",
    "    .squeeze()\n",
    ")\n",
    "ds = np.where(ds > 0, ds, np.nan)\n",
    "fig, ax = plt.subplots(figsize=(8, 8))\n",
    "\n",
    "ax.imshow(ds, cmap=\"tab20\")\n",
    "ax.set_axis_off()"
   ]
  }
 ],
 "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}