{ "cells": [ { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "## Proximity tools\n", "\n", "In this tutorial, you'll learn how to use [Xarray-Spatial's proximity toolset](https://xarray-spatial.readthedocs.io/en/latest/user_guide/proximity.html) to measure distances between points. You'll focus on an area of the Amazon rainforest around Oriximiná, Brazil, using data from the [European Commission Joint Research Centre's Global Surface Water Dataset](https://planetarycomputer.microsoft.com/dataset/jrc-gsw).\n", "\n", "Xarray-Spatial offers three proximity tools:\n", "* [Proximity Distance](https://xarray-spatial.readthedocs.io/en/latest/reference/_autosummary/xrspatial.proximity.proximity.html): for each point in the input raster, this tool calculates the distance to the nearest of a set of target points or source points.\n", "* [Proximity Allocation](https://xarray-spatial.readthedocs.io/en/latest/reference/_autosummary/xrspatial.proximity.allocation.html): for each cell in the input raster, this tool identifies the nearest source or target point (the 'allocation' point).\n", "* [Proximity Direction](https://xarray-spatial.readthedocs.io/en/latest/reference/_autosummary/xrspatial.proximity.direction.html): for each cell in the input raster, this tool returns the direction to the nearest source point (the 'allocation')." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import pystac_client\n", "import rasterio\n", "import rasterio.mask\n", "import numpy as np\n", "import planetary_computer\n", "import stackstac\n", "\n", "import matplotlib.pyplot as plt\n", "from matplotlib.colors import Normalize\n", "from matplotlib.colors import ListedColormap\n", "\n", "from xrspatial import proximity, allocation, direction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Prepare and load the data\n", "\n", "This example uses a small amount of data. Set up a local Dask \"cluster\" on a single machine to process the data in parallel." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2022-08-16 16:31:46,001 - distributed.diskutils - INFO - Found stale lock file and directory '/tmp/dask-worker-space/worker-4pshvxtk', purging\n", "2022-08-16 16:31:46,003 - distributed.diskutils - INFO - Found stale lock file and directory '/tmp/dask-worker-space/worker-bmfasp8e', purging\n", "2022-08-16 16:31:46,003 - distributed.diskutils - INFO - Found stale lock file and directory '/tmp/dask-worker-space/worker-ooa182jp', purging\n", "2022-08-16 16:31:46,003 - distributed.diskutils - INFO - Found stale lock file and directory '/tmp/dask-worker-space/worker-wey1_0cy', purging\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "/proxy/8787/status\n" ] } ], "source": [ "from dask.distributed import Client\n", "\n", "client = Client()\n", "print(f\"/proxy/{client.scheduler_info()['services']['dashboard']}/status\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You'll analyze a small area in the Amazon rainforest located around Oriximiná, State of Pará, Brazil. The region of interest contains parts of the Amazon River and smaller rivers such as the Nhamundá River, Trombetas River, and Paru de Oeste River. \n", "\n", "Use `pystac-client` to find all [STAC Items](https://github.com/radiantearth/stac-spec/blob/master/overview.md#item-overview) covering that area:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Returned 1 Items\n" ] } ], "source": [ "bounds = [-57.151965, -2.530125, -55.710724, -1.179033]\n", "\n", "catalog = pystac_client.Client.open(\n", " \"https://planetarycomputer.microsoft.com/api/stac/v1\",\n", " modifier=planetary_computer.sign_inplace,\n", ")\n", "jrc = catalog.search(collections=[\"jrc-gsw\"], bbox=bounds)\n", "\n", "items = list(jrc.item_collection())\n", "print(f\"Returned {len(items)} Items\")" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "item = items[0]\n", "item" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This notebook uses these three assets of the dataset: Maximum Water Extent, Water Seasonality, and Water Transitions. Use [stackstac.stack](https://stackstac.readthedocs.io/en/latest/api/main/stackstac.stack.html) to load the item data, select the desired assets, and crop the area to our bounds of interest. See the [jrc-gsw example notebook](../datasets/jrc-gsw/jrc-gsw-example.ipynb) to learn more about how to use this dataset." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'stackstac-e7dd0ef9e9571a9f4ddfcd90de32b185' (band: 3,\n",
       "                                                                y: 5405, x: 5766)>\n",
       "dask.array<getitem, shape=(3, 5405, 5766), dtype=float64, chunksize=(1, 3000, 3000), chunktype=numpy.ndarray>\n",
       "Coordinates: (12/16)\n",
       "    time            datetime64[ns] 2020-07-01\n",
       "    id              <U15 '60W_0Nv1_3_2020'\n",
       "  * band            (band) <U11 'extent' 'seasonality' 'transitions'\n",
       "  * x               (x) float64 -57.15 -57.15 -57.15 ... -55.71 -55.71 -55.71\n",
       "  * y               (y) float64 -1.179 -1.179 -1.179 -1.18 ... -2.53 -2.53 -2.53\n",
       "    sci:citation    <U169 'Jean-Francois Pekel, Andrew Cottam, Noel Gorelick,...\n",
       "    ...              ...\n",
       "    end_datetime    (band) object None '2020-12-31T11:59:59Z' None\n",
       "    proj:shape      object {40000}\n",
       "    proj:bbox       object {0.0, -60.0, -50.0, -10.0}\n",
       "    description     (band) <U189 'Provides information on all the locations e...\n",
       "    title           (band) <U20 'Maximum Water Extent' ... 'Water Transitions'\n",
       "    epsg            int64 4326\n",
       "Attributes:\n",
       "    spec:        RasterSpec(epsg=4326, bounds=(-57.152, -2.53025, -55.7105, -...\n",
       "    crs:         epsg:4326\n",
       "    transform:   | 0.00, 0.00,-57.15|\\n| 0.00,-0.00,-1.18|\\n| 0.00, 0.00, 1.00|\n",
       "    resolution:  0.00025
" ], "text/plain": [ "\n", "dask.array\n", "Coordinates: (12/16)\n", " time datetime64[ns] 2020-07-01\n", " id " ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "norm = Normalize(0, 255)\n", "\n", "fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(20, 15))\n", "\n", "for i, asset_key in enumerate(assets_of_interest):\n", " ax[i].imshow(\n", " data.sel(band=asset_key),\n", " norm=norm,\n", " cmap=cmaps[asset_key],\n", " )\n", " ax[i].set_title(asset_key, fontdict={\"fontsize\": 15})\n", " ax[i].set_axis_off()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Proximity Distance\n", "\n", "[xrspatial.proximity.proximity](https://xarray-spatial.readthedocs.io/en/latest/reference/_autosummary/xrspatial.proximity.proximity.html) takes a `values` raster and, for each point in the raster, computes the distance to the closest *target*. By default, all non-zero values in the raster are used as `target` pixels, but you can provide targets manually instead.\n", "\n", "The output raster is the same shape as the input, and the values distance to the closest target point (using the `distance_metric`, euclidean by default).\n", "\n", "By default, you'll use Euclidean as the `distance_metric`, the `target_values` parameter is all non-zero pixels in the values raster, and `max_distance` is set to infinity." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray '\\n    Water Extent proximity distance\\n    (Euclidean max_distance=infinity)\\n' (\n",
       "                                                                                                    y: 5405,\n",
       "                                                                                                    x: 5766)>\n",
       "array([[0.1372582 , 0.1370082 , 0.13675822, ..., 0.01077033, 0.01086566,\n",
       "        0.01096586],\n",
       "       [0.1372557 , 0.1370057 , 0.13675572, ..., 0.01053862, 0.01063602,\n",
       "        0.01073837],\n",
       "       [0.13725364, 0.13700365, 0.13675366, ..., 0.01030776, 0.01040733,\n",
       "        0.0105119 ],\n",
       "       ...,\n",
       "       [0.        , 0.        , 0.        , ..., 0.15853864, 0.15876181,\n",
       "        0.15898506],\n",
       "       [0.        , 0.        , 0.        , ..., 0.15865155, 0.15887456,\n",
       "        0.15909766],\n",
       "       [0.        , 0.        , 0.        , ..., 0.15876476, 0.15898761,\n",
       "        0.15921055]], dtype=float32)\n",
       "Coordinates: (12/16)\n",
       "    time            datetime64[ns] 2020-07-01\n",
       "    id              <U15 '60W_0Nv1_3_2020'\n",
       "    band            <U11 'extent'\n",
       "  * x               (x) float64 -57.15 -57.15 -57.15 ... -55.71 -55.71 -55.71\n",
       "  * y               (y) float64 -1.179 -1.179 -1.179 -1.18 ... -2.53 -2.53 -2.53\n",
       "    sci:citation    <U169 'Jean-Francois Pekel, Andrew Cottam, Noel Gorelick,...\n",
       "    ...              ...\n",
       "    end_datetime    object None\n",
       "    proj:shape      object {40000}\n",
       "    proj:bbox       object {0.0, -60.0, -50.0, -10.0}\n",
       "    description     <U189 'Provides information on all the locations ever det...\n",
       "    title           <U20 'Maximum Water Extent'\n",
       "    epsg            int64 4326\n",
       "Attributes:\n",
       "    spec:        RasterSpec(epsg=4326, bounds=(-57.152, -2.53025, -55.7105, -...\n",
       "    crs:         epsg:4326\n",
       "    transform:   | 0.00, 0.00,-57.15|\\n| 0.00,-0.00,-1.18|\\n| 0.00, 0.00, 1.00|\n",
       "    resolution:  0.00025
" ], "text/plain": [ "\n", "array([[0.1372582 , 0.1370082 , 0.13675822, ..., 0.01077033, 0.01086566,\n", " 0.01096586],\n", " [0.1372557 , 0.1370057 , 0.13675572, ..., 0.01053862, 0.01063602,\n", " 0.01073837],\n", " [0.13725364, 0.13700365, 0.13675366, ..., 0.01030776, 0.01040733,\n", " 0.0105119 ],\n", " ...,\n", " [0. , 0. , 0. , ..., 0.15853864, 0.15876181,\n", " 0.15898506],\n", " [0. , 0. , 0. , ..., 0.15865155, 0.15887456,\n", " 0.15909766],\n", " [0. , 0. , 0. , ..., 0.15876476, 0.15898761,\n", " 0.15921055]], dtype=float32)\n", "Coordinates: (12/16)\n", " time datetime64[ns] 2020-07-01\n", " id " ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(15, 7.5))\n", "\n", "extent_proximity_default.plot.imshow(cmap=\"Greens_r\", add_colorbar=False, ax=ax1)\n", "extent_data.plot.imshow(\n", " norm=norm, cmap=cmaps[\"extent\"], add_colorbar=False, alpha=0.5, ax=ax1\n", ")\n", "ax1.set_axis_off()\n", "ax1.set(title=extent_proximity_default.name)\n", "\n", "extent_proximity.plot.imshow(cmap=\"Greens_r\", add_colorbar=False, ax=ax2)\n", "extent_data.plot.imshow(\n", " norm=norm, cmap=cmaps[\"extent\"], add_colorbar=False, alpha=0.5, ax=ax2\n", ")\n", "ax2.set_axis_off()\n", "ax2.set(title=extent_proximity.name);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this visualization, dark green points are points that are closest to the surface water, while lighter points are further away. Blank points (all white) are further than `max_distance` to the closest point." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Proximity Allocation\n", "\n", "Similar to [xrspatial.proximity.proximity](https://xarray-spatial.readthedocs.io/en/latest/reference/_autosummary/xrspatial.proximity.proximity.html), [xrspatial.proximity.allocation](https://xarray-spatial.readthedocs.io/en/latest/reference/_autosummary/xrspatial.proximity.allocation.html) takes an aggregate as its input and finds the smallest distance from each cell to any one of the target points or source points. However, instead of returning the distance, [xrspatial.proximity.allocation](https://xarray-spatial.readthedocs.io/en/latest/reference/_autosummary/xrspatial.proximity.allocation.html) returns the value of the target point for each pixel.\n", "\n", "We'll compute the allocation, using `great_circle_distance` as our distance metric." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "3677.9271626991012" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from xrspatial import great_circle_distance\n", "\n", "d_great_circle = great_circle_distance(p1[0], p2[0], p1[1], p2[1])\n", "d_great_circle" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Compute the allocation proximity of the Water Transitions asset data:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "transitions_data = data.sel(band=\"transitions\").compute()\n", "\n", "transitions_allocation = allocation(\n", " transitions_data, distance_metric=\"GREAT_CIRCLE\", max_distance=d_great_circle\n", ")\n", "transitions_allocation.name = f\"\"\"\n", " Water Transitions proximity allocation\n", " (Great Circle max_distance={d_great_circle:.2f})\n", "\"\"\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Visualize the result:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[Text(0.5, 1.0, '\\n Water Transitions proximity allocation\\n (Great Circle max_distance=3677.93)\\n')]" ] }, "execution_count": 14, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(figsize=(10, 8))\n", "\n", "ax.set_axis_off()\n", "ax.imshow(transitions_allocation, norm=norm, cmap=cmaps[\"transitions\"])\n", "ax.imshow(transitions_data, norm=norm, cmap=cmaps[\"transitions\"], alpha=0.5)\n", "ax.set(title=transitions_allocation.name)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice the blocks in the visualization: each of the differently shaded blocks contains all of the points that share the target point in the center as their nearest target point." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Proximity Direction\n", "\n", "Finally, we'll use [xrspatial.proximity.direction](https://xarray-spatial.readthedocs.io/en/latest/reference/_autosummary/xrspatial.proximity.direction.html) to compute the *direction* to the nearest target, rather than the distance. The output values range from 0 to 360, with\n", "\n", "- 0 is the source cell itself\n", "- 90 is East\n", "- 180 is South\n", "- 270 is West\n", "- 360 is North\n", "\n", "In this example, we'll use the `seasonality` band, the `manhattan_distance` metric, and set `target_values=[10, 11, 12]`. The values in `seasonality` indicate the number of months per year that surface water is present in a pixel, so we'll find direction to the nearest pixel with water coverage for at least 10 months of the year." ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.04169199999999895" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from xrspatial import manhattan_distance\n", "\n", "d_manhattan = manhattan_distance(p1[0], p2[0], p1[1], p2[1])\n", "d_manhattan" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Calculate the proximity direction for the Water Seasonality asset data:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "seasonality_direction = direction(\n", " data.sel(band=\"seasonality\"),\n", " distance_metric=\"MANHATTAN\",\n", " target_values=[10, 11, 12],\n", " max_distance=d_manhattan,\n", ").compute()\n", "\n", "seasonality_direction.name = f\"\"\"\n", " 10-12 months Water Seasonality proximity direction\n", " (Manhattan max_distance={d_manhattan:.5f})\n", "\"\"\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To better visualization all targets, filter out non-target pixels:" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [], "source": [ "seasonality_data = (\n", " data.sel(band=\"seasonality\")\n", " .where(lambda x: x >= 10, other=0)\n", " .rename(\"Water Seasonality (10-12 months)\")\n", ").compute()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Plot the seasonality data and the direction proximity data next to each other." ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "data": { "text/html": [ "" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, (ax1, ax2) = plt.subplots(ncols=2, figsize=(15, 7.5))\n", "\n", "ax1.set_axis_off()\n", "ax2.set_axis_off()\n", "\n", "seasonality_data.plot.imshow(\n", " norm=norm, cmap=cmaps[\"seasonality\"], ax=ax1, add_colorbar=False\n", ")\n", "ax1.set(title=seasonality_data.name)\n", "\n", "seasonality_direction.plot.imshow(cmap=\"Blues\", ax=ax2, add_colorbar=False)\n", "seasonality_data.plot.imshow(\n", " norm=norm, cmap=cmaps[\"seasonality\"], alpha=0.5, ax=ax2, add_colorbar=False\n", ")\n", "ax2.set(title=seasonality_direction.name);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Next steps: analyze more datasets\n", "\n", "Try using Xarray-Spatial's [Proximity Distance](https://xarray-spatial.readthedocs.io/en/latest/user_guide/proximity.html#Proximity-Distance), [Proximity Allocation](https://xarray-spatial.readthedocs.io/en/latest/user_guide/proximity.html#Proximity-Allocation), and [Proximity Direction](https://xarray-spatial.readthedocs.io/en/latest/user_guide/proximity.html#Proximity-Direction) functions with different areas of interest and different assets of the [European Commission Joint Research Centre's Global Surface Water Dataset](https://planetarycomputer.microsoft.com/dataset/jrc-gsw). Or try using these tools with the [USGS GAP/LANDFIRE National Terrestrial Ecosystems dataset](https://planetarycomputer.microsoft.com/dataset/gap) to analyze proximity distances between different kinds of landcovers." ] } ], "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": { "cc4ca06f56dc49ebbce15e9dd6c2696a": { "model_module": "@jupyter-widgets/base", "model_module_version": "1.2.0", "model_name": "LayoutModel", "state": {} }, "fc2f6a422d2144a385ff318fffeefe63": { "model_module": "@jupyter-widgets/controls", "model_module_version": "1.5.0", "model_name": "VBoxModel", "state": { "layout": "IPY_MODEL_cc4ca06f56dc49ebbce15e9dd6c2696a" } } }, "version_major": 2, "version_minor": 0 } } }, "nbformat": 4, "nbformat_minor": 4 }