{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Focal Statistics\n", "\n", "In this tutorial we calculate focal statistics and determine hot and cold spots of NDVI (Normalized difference vegetation index) over a time series of satellite images. On its own, NDVI is used to highlight live green vegetation. Its hot and cold spots help determine the growth or loss of plants. In this notebook, we'll see how to:\n", "\n", "- Search for satellite data by item ID using `pystac_client`\n", "- Visualize true color images\n", "- Calculate NDVI\n", "- Smooth images with a mean filter\n", "- Calculate focal statistics of the values within a specified focal neighborhood for each pixel in an input data array\n", "- Identify hot and cold spots in a image, neighborhoods that are significantly different from the rest of the image\n", "\n", "The focus of this notebook is to analyse information for each pixel based on its focal neighborhood kernel. The [xrspatial.focal](https://xarray-spatial.readthedocs.io/en/latest/user_guide/focal.html) module from `xarray-spatial` provides a set of analysis tools performing neighborhood operations that will be used through this tutorial." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import xarray as xr\n", "\n", "import stackstac\n", "import planetary_computer\n", "import pystac_client\n", "\n", "import matplotlib.pyplot as plt\n", "\n", "import xrspatial.multispectral as ms\n", "from xrspatial.convolution import calc_cellsize, circle_kernel, convolution_2d\n", "from xrspatial.focal import mean, focal_stats, hotspots\n", "\n", "from dask.distributed import Client, progress" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Local Dask Cluster\n", "\n", "We'll use a small number of images for this example. We'll parallelize reading the data from Azure Blob Storage using a local Dask \"cluster\" on this single machine." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "tags": [] }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "2022-08-16 16:30:24,531 - distributed.diskutils - INFO - Found stale lock file and directory '/tmp/dask-worker-space/worker-l6lzssem', purging\n", "2022-08-16 16:30:24,532 - distributed.diskutils - INFO - Found stale lock file and directory '/tmp/dask-worker-space/worker-h374kdv4', purging\n", "2022-08-16 16:30:24,532 - distributed.diskutils - INFO - Found stale lock file and directory '/tmp/dask-worker-space/worker-isifihn_', purging\n", "2022-08-16 16:30:24,532 - distributed.diskutils - INFO - Found stale lock file and directory '/tmp/dask-worker-space/worker-swbsk4a4', purging\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "/proxy/8787/status\n" ] } ], "source": [ "client = Client()\n", "print(f\"/proxy/{client.scheduler_info()['services']['dashboard']}/status\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can access the Dask Dashboard by pasting that URL into the Dask labextension field. See [Scale with Dask](../quickstarts/scale-with-dask.ipynb) for more." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The region of interest is a small area in the Amazon rainforest located in State of Mato Grosso and State of Amazonas, Brazil. In order to calculate NDVI accurately, we found these least cloudy scenes by searching with the [STAC API](https://planetarycomputer.microsoft.com/docs/quickstarts/reading-stac/) and filtering the results." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "catalog = pystac_client.Client.open(\n", " \"https://planetarycomputer.microsoft.com/api/stac/v1\",\n", " modifier=planetary_computer.sign_inplace,\n", ")\n", "\n", "ids = [\n", " \"S2A_MSIL2A_20200616T141741_R010_T20LQR_20200822T232052\",\n", " \"S2B_MSIL2A_20190617T142049_R010_T20LQR_20201006T032921\",\n", " \"S2B_MSIL2A_20180712T142039_R010_T20LQR_20201011T150557\",\n", " \"S2B_MSIL2A_20170727T142039_R010_T20LQR_20210210T153028\",\n", " \"S2A_MSIL2A_20160627T142042_R010_T20LQR_20210211T234456\",\n", "]\n", "search = catalog.search(collections=[\"sentinel-2-l2a\"], ids=ids)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we'll sign the STAC items so we can download the data from blob storage. See [Using Tokens for Data Access](../concepts/sas.ipynb) for more. " ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
<xarray.DataArray 'stackstac-6aab09147bacb75bf1d5e6bdb3733648' (time: 5,\n", " band: 4,\n", " y: 227, x: 226)>\n", "dask.array<where, shape=(5, 4, 227, 226), dtype=float64, chunksize=(1, 1, 227, 226), chunktype=numpy.ndarray>\n", "Coordinates: (12/46)\n", " * time (time) datetime64[ns] 2016-06-28...\n", " id (time) <U54 'S2A_MSIL2A_20160627...\n", " * band (band) <U5 'blue' 'green' ... 'nir'\n", " * x (x) float64 1.362e+06 ... 1.474e+06\n", " * y (y) float64 -9.075e+05 ... -1.02...\n", " s2:granule_id (time) <U62 'S2A_OPER_MSI_L2A_TL...\n", " ... ...\n", " title (band) <U20 'Band 2 - Blue - 10m...\n", " proj:transform object {0.0, 9100000.0, 10.0, -1...\n", " common_name (band) <U5 'blue' 'green' ... 'nir'\n", " center_wavelength (band) float64 0.49 0.56 ... 0.842\n", " full_width_half_max (band) float64 0.098 ... 0.145\n", " epsg int64 32619\n", "Attributes:\n", " spec: RasterSpec(epsg=32619, bounds=(1361500, -1021000, 1474500, -...\n", " crs: epsg:32619\n", " transform: | 500.00, 0.00, 1361500.00|\\n| 0.00,-500.00,-907500.00|\\n| 0...\n", " resolution: 500