{ "cells": [ { "cell_type": "markdown", "id": "4d70b9e3-6d9f-4eb0-961d-3d2d8889541a", "metadata": {}, "source": [ "## Using ALOS PALSAR and Forest/Non-Forest Annual Mosaics with the Planetary Computer STAC API" ] }, { "cell_type": "markdown", "id": "c8ec9787-3870-47d3-8b0c-38e2407e08b9", "metadata": {}, "source": [ "ALOS PALSAR Annual Mosaic consists of radar backscatter in the HH and HV with additional bands for ancillary information about incidince angle, masking, and date of acquisition (per pixel).\n", "\n", "For this example we'll find a 1x1 degree data tile with a diversity of values, and plot the HH, HV and matching Forest Classification." ] }, { "cell_type": "code", "execution_count": 1, "id": "051fa411-8bdb-43aa-8572-1b17262a1b70", "metadata": { "tags": [] }, "outputs": [], "source": [ "import numpy as np\n", "import planetary_computer\n", "import stackstac\n", "import pystac_client\n", "import warnings\n", "\n", "warnings.filterwarnings(\"ignore\", \"The argument 'infer_datetime_format'\", UserWarning)" ] }, { "cell_type": "markdown", "id": "559bb25e-b9ba-4c1f-b77b-8044268a635d", "metadata": {}, "source": [ "### Polarization Returns\n", "\n", "We'll use the Planetary Computer's STAC API to find some scenes matching our area of interest. See [reading STAC](https://planetarycomputer.microsoft.com/docs/quickstarts/reading-stac/) for more." ] }, { "cell_type": "code", "execution_count": 2, "id": "2607dea7-7698-47b2-a4fc-6acbc52ec41f", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Returned 21 items\n" ] } ], "source": [ "bbox = [9.4, 0, 9.5, 1]\n", "client = pystac_client.Client.open(\n", " \"https://planetarycomputer.microsoft.com/api/stac/v1/\",\n", " modifier=planetary_computer.sign_inplace,\n", ")\n", "search = client.search(\n", " collections=[\"alos-palsar-mosaic\"],\n", " bbox=bbox,\n", ")\n", "items = search.item_collection()\n", "print(f\"Returned {len(items)} items\")" ] }, { "cell_type": "markdown", "id": "e6657b7d-0992-45f3-919e-c79c662092b8", "metadata": {}, "source": [ "Each of these items has a handful of assets." ] }, { "cell_type": "code", "execution_count": 3, "id": "7f205ba6-a343-4d6d-a7f2-c37afa43547b", "metadata": { "tags": [] }, "outputs": [ { "data": { "text/plain": [ "['HH',\n", " 'HV',\n", " 'date',\n", " 'mask',\n", " 'linci',\n", " 'metadata',\n", " 'tilejson',\n", " 'rendered_preview']" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "item = items[0]\n", "list(item.assets)" ] }, { "cell_type": "markdown", "id": "859b7032-e5ff-4465-b0c7-3fd070a1d0bc", "metadata": {}, "source": [ "We'll select the first item and load it into xarray using stackstac. Note that the underlying data have a `uint16` dtype, so we specify that when creating the `DataArray`." ] }, { "cell_type": "code", "execution_count": 4, "id": "bd018379-521b-4ebc-8431-0d0819f3a50a", "metadata": { "tags": [] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'stackstac-1a4421f08bc5ca33d576ce09d60ce330' (band: 2,\n",
       "                                                                y: 4501, x: 4500)>\n",
       "dask.array<getitem, shape=(2, 4501, 4500), dtype=uint16, chunksize=(1, 1024, 1024), chunktype=numpy.ndarray>\n",
       "Coordinates: (12/27)\n",
       "    time                            datetime64[ns] 2021-01-01\n",
       "    id                              <U21 'N02E009_21_F02DAR_MOS'\n",
       "  * band                            (band) <U2 'HH' 'HV'\n",
       "  * x                               (x) float64 9.0 9.0 9.0 ... 9.999 10.0 10.0\n",
       "  * y                               (y) float64 2.0 2.0 2.0 ... 1.001 1.0 1.0\n",
       "    proj:epsg                       int64 4326\n",
       "    ...                              ...\n",
       "    start_datetime                  <U20 '2021-01-01T00:00:00Z'\n",
       "    platform                        <U6 'ALOS-2'\n",
       "    gsd                             int64 25\n",
       "    proj:bbox                       object {9.0, 10.0, 2.0000000000000004, 1....\n",
       "    raster:bands                    object {'nodata': 1, 'data_type': 'uint16'}\n",
       "    epsg                            int64 4326\n",
       "Attributes:\n",
       "    spec:        RasterSpec(epsg=4326, bounds=(9.0, 1.0, 10.0, 2.000222222222...\n",
       "    crs:         epsg:4326\n",
       "    transform:   | 0.00, 0.00, 9.00|\\n| 0.00,-0.00, 2.00|\\n| 0.00, 0.00, 1.00|\n",
       "    resolution:  0.00022222222222222223
" ], "text/plain": [ "\n", "dask.array\n", "Coordinates: (12/27)\n", " time datetime64[ns] 2021-01-01\n", " id " ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "data.sel(band=\"HH\").plot.hist(bins=100, size=5);" ] }, { "cell_type": "code", "execution_count": 6, "id": "698b7bf4-65f1-41a1-9777-f96a249cbe1a", "metadata": { "tags": [] }, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "img = data.sel(band=\"HV\").plot(vmin=0, vmax=14000, size=8)\n", "img.axes.set_axis_off();" ] }, { "cell_type": "markdown", "id": "ba6be223-8edf-4592-8b45-423470ecf5d1", "metadata": {}, "source": [ "### Forest Classification\n", "\n", "ALOS Forest/Non-Forest Classification is derived from the ALOS PALSAR Annual Mosaic, and classifies the pixels to detect forest cover. We'll search for Forest / Non-Forest items covering the same area." ] }, { "cell_type": "code", "execution_count": 7, "id": "f36ce133-90c7-4fe0-a472-9517fae4d8a8", "metadata": { "tags": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Returned 6 items\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/srv/conda/envs/notebook/lib/python3.10/site-packages/pystac_client/item_search.py:841: FutureWarning: get_all_items() is deprecated, use item_collection() instead.\n", " warnings.warn(\n" ] } ], "source": [ "import shapely.geometry\n", "\n", "s = shapely.geometry.shape(item.geometry)\n", "intersects = shapely.geometry.mapping(s.centroid)\n", "\n", "search = client.search(collections=[\"alos-fnf-mosaic\"], intersects=intersects)\n", "fnf_items = search.get_all_items()\n", "print(f\"Returned {len(fnf_items)} items\")" ] }, { "cell_type": "markdown", "id": "fbf91d1c-92fe-4ed7-bd58-7d80f1264950", "metadata": {}, "source": [ "The primary asset is under the key `C`." ] }, { "cell_type": "code", "execution_count": 8, "id": "46bc3703-510a-41ab-be41-bf3642313a8d", "metadata": { "tags": [] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'stackstac-3882cecac9b4ee67fd63e1c7d3c237f3' (y: 4500, x: 4500)>\n",
       "dask.array<getitem, shape=(4500, 4500), dtype=uint8, chunksize=(1024, 1024), chunktype=numpy.ndarray>\n",
       "Coordinates: (12/18)\n",
       "    time            datetime64[ns] 2020-01-01\n",
       "    id              <U14 'N02E009_20_FNF'\n",
       "    band            <U1 'C'\n",
       "  * x               (x) float64 9.0 9.0 9.0 9.001 ... 9.999 9.999 10.0 10.0\n",
       "  * y               (y) float64 2.0 2.0 2.0 1.999 1.999 ... 1.001 1.001 1.0 1.0\n",
       "    proj:shape      object {4500}\n",
       "    ...              ...\n",
       "    instruments     <U8 'PALSAR-2'\n",
       "    title           <U1 'C'\n",
       "    gsd             int64 25\n",
       "    proj:bbox       object {9.0, 10.0, 2.0, 1.0}\n",
       "    raster:bands    object {'nodata': 0, 'data_type': 'uint8'}\n",
       "    epsg            int64 4326\n",
       "Attributes:\n",
       "    spec:        RasterSpec(epsg=4326, bounds=(9.0, 1.0, 10.0, 2.0), resoluti...\n",
       "    crs:         epsg:4326\n",
       "    transform:   | 0.00, 0.00, 9.00|\\n| 0.00,-0.00, 2.00|\\n| 0.00, 0.00, 1.00|\n",
       "    resolution:  0.00022222222222222223
" ], "text/plain": [ "\n", "dask.array\n", "Coordinates: (12/18)\n", " time datetime64[ns] 2020-01-01\n", " id " ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import matplotlib.colors\n", "\n", "cmap = matplotlib.colors.ListedColormap(\n", " [\n", " (0, 0, 0, 0), # nodata\n", " (0, 0.7, 0, 1), # forest, > 90% canopy\n", " (0.51, 0.94, 0.38, 1), # forest 10-90% canopy\n", " (1, 1, 0.6, 1), # non-forest\n", " (0.0, 0.0, 1, 1), # water\n", " ]\n", ")\n", "\n", "# Plot the classification map\n", "p = fnf_tile.plot(cmap=cmap, vmin=0, vmax=4, size=8)\n", "ticks = np.arange(0.4, 4, 0.8)\n", "p.colorbar.set_ticks(ticks, labels=labels)\n", "p.colorbar.set_label(\"Class\")" ] } ], "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.10" }, "widgets": { "application/vnd.jupyter.widget-state+json": { "state": { "a8c4d99d835b44bfb4cc62c89f88abc8": { "model_module": "@jupyter-widgets/base", "model_module_version": "1.2.0", "model_name": "LayoutModel", "state": {} }, "c5649c33bf18455e88ab68b9c8e5e8b3": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "VBoxModel", "state": { "layout": "IPY_MODEL_a8c4d99d835b44bfb4cc62c89f88abc8" } } }, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 5 }