{ "cells": [ { "cell_type": "markdown", "id": "d87f2c82", "metadata": {}, "source": [ "## Accessing NAIP data with the Planetary Computer STAC API\n", "\n", "The [National Agriculture Imagery Program](https://www.fsa.usda.gov/programs-and-services/aerial-photography/imagery-programs/naip-imagery/) (NAIP) provides US-wide, high-resolution aerial imagery, with four spectral bands (R, G, B, IR). NAIP is administered by the [Aerial Field Photography Office](https://www.fsa.usda.gov/programs-and-services/aerial-photography/) (AFPO) within the [US Department of Agriculture](https://www.usda.gov/) (USDA). Data are captured at least once every three years for each state. This dataset represents NAIP data from 2010-present, in [cloud-optimized GeoTIFF](https://www.cogeo.org/) format.\n", "\n", "This notebook demonstrates the use of the Planetary Computer STAC API to query for NAIP imagery." ] }, { "cell_type": "markdown", "id": "e86824ea", "metadata": {}, "source": [ "### 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": "0062a353", "metadata": {}, "outputs": [], "source": [ "import pystac_client\n", "import planetary_computer\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()" ] }, { "cell_type": "markdown", "id": "97edd22e-7886-4734-9fb5-7ec4038cd164", "metadata": {}, "source": [ "### Data access\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, "id": "2a343e0d-9931-4f38-922b-d2d8e8cc8645", "metadata": {}, "outputs": [], "source": [ "catalog = pystac_client.Client.open(\n", " \"https://planetarycomputer.microsoft.com/api/stac/v1\",\n", " modifier=planetary_computer.sign_inplace,\n", ")" ] }, { "cell_type": "markdown", "id": "47354901", "metadata": {}, "source": [ "### Choose our region and times of interest\n", "\n", "This area is in South Jordan, Utah, which the Internet says is one of the fastest-growing cities in the US. Let's see whether we can see some development in this area in the time spanned by our NAIP collection." ] }, { "cell_type": "code", "execution_count": 3, "id": "6deaf8e8", "metadata": {}, "outputs": [], "source": [ "area_of_interest = {\n", " \"type\": \"Polygon\",\n", " \"coordinates\": [\n", " [\n", " [-111.9839859008789, 40.5389819819361],\n", " [-111.90502166748045, 40.5389819819361],\n", " [-111.90502166748045, 40.57015381856105],\n", " [-111.9839859008789, 40.57015381856105],\n", " [-111.9839859008789, 40.5389819819361],\n", " ]\n", " ],\n", "}" ] }, { "cell_type": "markdown", "id": "3a38a21c", "metadata": {}, "source": [ "This collection includes data from 2010, so we'll search for one image near the beginning of that range and one from 2018." ] }, { "cell_type": "code", "execution_count": 4, "id": "051121fa", "metadata": {}, "outputs": [], "source": [ "range_old = \"2010-01-01/2013-01-01\"\n", "range_new = \"2018-01-01/2021-01-01\"" ] }, { "cell_type": "code", "execution_count": 5, "id": "584b2976", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "4 Items found in the 'old' range\n", "4 Items found in the 'new' range\n" ] } ], "source": [ "search_old = catalog.search(\n", " collections=[\"naip\"], intersects=area_of_interest, datetime=range_old\n", ")\n", "\n", "search_new = catalog.search(\n", " collections=[\"naip\"], intersects=area_of_interest, datetime=range_new\n", ")\n", "\n", "items_old = search_old.item_collection()\n", "items_new = search_new.item_collection()\n", "\n", "print(f\"{len(items_old)} Items found in the 'old' range\")\n", "print(f\"{len(items_new)} Items found in the 'new' range\")" ] }, { "cell_type": "markdown", "id": "fcb22cc7", "metadata": {}, "source": [ "As seen above, there are multiple items that intersect our area of interest for each year. The following code will choose the item that has the most overlap:" ] }, { "cell_type": "code", "execution_count": 6, "id": "a9bd9f03", "metadata": {}, "outputs": [], "source": [ "from shapely.geometry import shape\n", "\n", "area_shape = shape(area_of_interest)\n", "target_area = area_shape.area\n", "\n", "\n", "def area_of_overlap(item):\n", " overlap_area = shape(item.geometry).intersection(shape(area_of_interest)).area\n", " return overlap_area / target_area\n", "\n", "\n", "item_old = sorted(items_old, key=area_of_overlap, reverse=True)[0]\n", "item_new = sorted(items_new, key=area_of_overlap, reverse=True)[0]" ] }, { "cell_type": "markdown", "id": "1bab7193", "metadata": {}, "source": [ "### Render images\n", "\n", "Each Item has a `rendered_preview` which uses the Planetary Computer's Data API to dynamically create a preview image from the raw data." ] }, { "cell_type": "code", "execution_count": 7, "id": "8c46a5a4-bdbf-401d-bb4d-f89146a19fa5", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.display import Image\n", "\n", "Image(url=item_old.assets[\"rendered_preview\"].href)" ] }, { "cell_type": "code", "execution_count": 8, "id": "a8864522-0c93-42f8-bc22-2af72b29555b", "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "Image(url=item_new.assets[\"rendered_preview\"].href)" ] }, { "cell_type": "markdown", "id": "22ec75f7-ec8c-4660-8421-d88ec14800c6", "metadata": {}, "source": [ "To read the raw Cloud Optimized GeoTIFF, use a library like [rioxarray](https://corteva.github.io/rioxarray/html/rioxarray.html) or [rasterio](https://rasterio.readthedocs.io/) and the `image` asset." ] }, { "cell_type": "code", "execution_count": 9, "id": "0ab9b2d6", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (band: 3, y: 7630, x: 6000)>\n",
       "[137340000 values with dtype=uint8]\n",
       "Coordinates:\n",
       "  * band         (band) int64 1 2 3\n",
       "  * x            (x) float64 4.15e+05 4.15e+05 4.15e+05 ... 4.209e+05 4.209e+05\n",
       "  * y            (y) float64 4.491e+06 4.491e+06 ... 4.483e+06 4.483e+06\n",
       "    spatial_ref  int64 0\n",
       "Attributes:\n",
       "    AREA_OR_POINT:  Area\n",
       "    scale_factor:   1.0\n",
       "    add_offset:     0.0
" ], "text/plain": [ "\n", "[137340000 values with dtype=uint8]\n", "Coordinates:\n", " * band (band) int64 1 2 3\n", " * x (x) float64 4.15e+05 4.15e+05 4.15e+05 ... 4.209e+05 4.209e+05\n", " * y (y) float64 4.491e+06 4.491e+06 ... 4.483e+06 4.483e+06\n", " spatial_ref int64 0\n", "Attributes:\n", " AREA_OR_POINT: Area\n", " scale_factor: 1.0\n", " add_offset: 0.0" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import rioxarray\n", "\n", "ds = rioxarray.open_rasterio(item_old.assets[\"image\"].href).sel(band=[1, 2, 3])\n", "ds" ] } ], "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.10.6" }, "widgets": { "application/vnd.jupyter.widget-state+json": { "state": {}, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 5 }