datasets/3dep/3dep-seamless-example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Accessing 3DEP Seamless DEM data with the Planetary Computer STAC API"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Environment setup\n",
    "\n",
    "This notebook works with or without a Planetary Computer subscription key, but you will be given more permissive access to the data with a key. The Planetary Computer Hub is pre-configured to use your key."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import planetary_computer\n",
    "import pystac_client"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Open and explore the 3dep-seamless collection\n",
    "\n",
    "[3DEP](https://planetarycomputer.microsoft.com/dataset/3dep-seamless) is a national elevation dataset program from the [USGS](https://www.usgs.gov/core-science-systems/ngp/3dep).\n",
    "It provides elevation data for the entire United States using a wide variety of survey methods, including modern lidar airborne surveys.\n",
    "The Planetary Computer hosts seamless elevation Digital Elevation Models (DEMs) at two resolutions, 1/3 arc-second and 1 arc-second.\n",
    "\n",
    "The datasets hosted by the Planetary Computer are available from [Azure Blob Storage](https://docs.microsoft.com/en-us/azure/storage/blobs/). We'll use [pystac-client](https://pystac-client.readthedocs.io/) to search the Planetary Computer's [STAC API](https://planetarycomputer.microsoft.com/api/stac/v1/docs) for the subset of the data that we care about, and then we'll load the data directly from Azure Blob Storage. We'll specify a `modifier` so that we can access the data stored in the Planetary Computer's private Blob Storage Containers. See [Reading from the STAC API](https://planetarycomputer.microsoft.com/docs/quickstarts/reading-stac/) and [Using tokens for data access](https://planetarycomputer.microsoft.com/docs/concepts/sas/) for more."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "catalog = pystac_client.Client.open(\n",
    "    \"https://planetarycomputer.microsoft.com/api/stac/v1\",\n",
    "    modifier=planetary_computer.sign_inplace,\n",
    ")\n",
    "threedep = catalog.get_child(id=\"3dep-seamless\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The 3DEP seamless dataset covers the entire United States, as you can see from the bounds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'bbox': [[-174.001666666983,\n",
       "   -15.00166666667,\n",
       "   164.0016666666,\n",
       "   84.00166666666]]}"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "threedep.extent.spatial.to_dict()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plotting elevation data\n",
    "\n",
    "Let's take a look at Mt. Shasta in Northern California.\n",
    "We define a point, then find a bounding box around that point.\n",
    "The query area for this example is pretty small, and we will only fetch one tile — you can increase the buffer size and this example should work with multiple tiles (though the plotting will take longer)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2 items found\n"
     ]
    }
   ],
   "source": [
    "latitude = 41.41505315463394\n",
    "longitude = -122.19468933831948\n",
    "buffer = 0.1\n",
    "bbox = [longitude - buffer, latitude - buffer, longitude + buffer, latitude + buffer]\n",
    "search = catalog.search(collections=\"3dep-seamless\", bbox=bbox)\n",
    "items = list(search.item_collection())\n",
    "print(f\"{len(items)} items found\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Planetary Computer hosts both the 1 arc-second (nominal 30m pixel size) and 1/3 arc-second (nominal 10m pixel size) resolution 3DEP data.\n",
    "Let's separate the two sets into their own lists by filtering on the Ground Sample Distance (GSD) in the STAC items.\n",
    "\n",
    "We'll also [sign the assets](../quickstarts/reading-stac.ipynb) before downloading, which can be done with or without a Planetary Computer subscription key."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "items_low_res = [item.to_dict() for item in items if item.properties[\"gsd\"] == 30]\n",
    "items_high_res = [item.to_dict() for item in items if item.properties[\"gsd\"] == 10]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "DEMs are relatively simple raster datasets — they only have one band (elevation). Let's compare visualizations of high and low resolution images for our area of interest.\n",
    "We'll read the items into a DatayArray using [`stackstac`](https://stackstac.readthedocs.io/), taking care to crop the larger assets down to our area of interest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "import stackstac\n",
    "import warnings\n",
    "\n",
    "warnings.filterwarnings(\"ignore\", message=\"The argument\")\n",
    "\n",
    "low_res_data = stackstac.stack(items_low_res, bounds=bbox).squeeze().compute()\n",
    "high_res_data = stackstac.stack(items_high_res, bounds=bbox).squeeze().compute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"https://ai4edatasetspublicassets.blob.core.windows.net/assets/notebook-output/.-datasets-3dep-3dep-seamless-example.ipynb/14.png\"/>"
      ],
      "text/plain": [
       "<Figure size 2400x1200 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, axes = plt.subplots(ncols=2, figsize=(24, 12), sharex=True, sharey=True)\n",
    "low_res_data.plot.imshow(ax=axes[0], add_colorbar=False)\n",
    "high_res_data.plot.imshow(ax=axes[1], add_colorbar=False)\n",
    "\n",
    "axes[0].set(title=\"1 arc-second\")\n",
    "axes[1].set(title=r\"$\\frac{1}{3}$ arc-second\");"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "a0bc0066fb6b219d769b3dc267816f1eb584a63f4b10cfd1c205fe8df8fe8802"
  },
  "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.11.8"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {
     "550e9c6c01c647a9afe71d92d7471681": {
      "model_module": "@jupyter-widgets/base",
      "model_module_version": "2.0.0",
      "model_name": "LayoutModel",
      "state": {}
     },
     "821e2fa4ecba410ab745f3f401f362c4": {
      "model_module": "@jupyter-widgets/controls",
      "model_module_version": "2.0.0",
      "model_name": "VBoxModel",
      "state": {
       "layout": "IPY_MODEL_550e9c6c01c647a9afe71d92d7471681"
      }
     }
    },
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}