{ "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()\n", "\n", "client = pystac_client.Client.open(\n", " \"https://planetarycomputer.microsoft.com/api/stac/v1/\",\n", " modifier=planetary_computer.sign_inplace,\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': ,\n", " '3dep-lidar-pointsourceid': ,\n", " '3dep-lidar-intensity': ,\n", " '3dep-lidar-hag': ,\n", " '3dep-lidar-dtm-native': ,\n", " '3dep-lidar-dtm': ,\n", " '3dep-lidar-dsm': ,\n", " '3dep-lidar-classification': }" ] }, "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 = {x.collection_id: x for x in search.get_all_items()}\n", "items" ] }, { "cell_type": "code", "execution_count": 3, "id": "1761c840-031a-41ed-b74e-0b673a7fa2ce", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'https://usgslidareuwest.blob.core.windows.net/usgs-3dep-cogs/usgs-cogs/USGS_LPC_VA_Fairfax_County_2018/numberofreturns/USGS_LPC_VA_Fairfax_County_2018-returns-5m-2-1.tif?st=2022-08-08T14%3A32%3A43Z&se=2022-08-16T14%3A32%3A43Z&sp=rl&sv=2021-06-08&sr=c&skoid=c85c15d6-d1ae-42d4-af60-e2ca0f81359b&sktid=72f988bf-86f1-41af-91ab-2d7cd011db47&skt=2022-08-09T14%3A32%3A42Z&ske=2022-08-16T14%3A32%3A42Z&sks=b&skv=2021-06-08&sig=pNxIEF9eBA0rbTbQIV5i/XqtwzBzzqMHOw7J7JWR3TA%3D'" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "items[\"3dep-lidar-returns\"].assets[\"data\"].href" ] }, { "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": 4, "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": 5, "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": 6, "id": "29c880bf-b0e8-4be7-9622-98c4726e98b3", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "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": 7, "id": "b6cd8021-4257-46e4-98fe-9f286445fe46", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "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": 8, "id": "c6825644-d141-4204-baa3-477366c45dae", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "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": 9, "id": "a154eafe-dc47-472a-98f2-78ff0ee7563c", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "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": 10, "id": "4f7484c0-8ec7-4428-ae76-6fb91c13180c", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "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": 11, "id": "3afbef14-f6e7-4ab5-a7aa-20b2801c1914", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "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": 12, "id": "48919fd8-39d0-4edd-ae65-61ce829771b7", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "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": 13, "id": "ae507018-d725-482c-b90e-1cebd5bc7741", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "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.13" }, "widgets": { "application/vnd.jupyter.widget-state+json": { "state": {}, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 5 }