{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Examples\n", "\n", "We'll watch coastline movement from [longshore drift](https://en.wikipedia.org/wiki/Longshore_drift) on [Cape Cod](https://www.google.com/maps/place/Chatham,+MA/@41.7498076,-70.2026227,10.73z/data=!4m13!1m7!3m6!1s0x89fb15440149e94d:0x1f9c0efa001cb20b!2sCape+Cod!3b1!8m2!3d41.6687897!4d-70.2962408!3m4!1s0x89fb142168afbe53:0x714436ec7d485a53!8m2!3d41.6821432!4d-69.9597359), Massachusetts, USA.\n", "\n", "We'll make quarterly, cloud-free-ish composites of Sentinel-2 imagery from AWS. Then turn them into GIFs." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import pystac_client\n", "import stackstac\n", "import dask.array as da\n", "import distributed\n", "from geogif import gif, dgif" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Launch a local Dask cluster. (This would be much faster with a cluster in the cloud, FWIW.)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ ], "source": [ "client = distributed.Client()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Search for STAC items overlapping our area of interest, from AWS's open data catalog." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "bbox = (-70.05487089045347, 41.59659777453696, -69.8481905926019, 41.69921590764708)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 1.41 s, sys: 167 ms, total: 1.57 s\n", "Wall time: 5.9 s\n" ] } ], "source": [ "%%time\n", "items = (\n", " pystac_client.Client.open(\"https://earth-search.aws.element84.com/v1\")\n", " .search(\n", " bbox=bbox,\n", " collections=[\"sentinel-2-l2a\"],\n", " datetime=\"2018-01-01/2024-01-01\",\n", " )\n", ").item_collection()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Make quarterly median composites\n", "\n", "Turn STAC items into xarray as a temporal stack, using [stackstac](https://stackstac.readthedocs.io/).\n", "\n", "Then mask out bad (cloudy) pixels, according to the Sentinel-2 SCL Scene Classification Map, and take the temporal median of each quarter (three months) to hopefully get an okay-looking cloud-free frame representative of those three months.\n", "\n", "Each 3-month composite will be one frame in our animation. We end up with 18 frames." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
<xarray.DataArray 'stackstac-706a4f359b3be98237baef766ab69231' (time: 24,\n",
" band: 3,\n",
" y: 387, x: 579)>\n",
"dask.array<stack, shape=(24, 3, 387, 579), dtype=float64, chunksize=(1, 1, 387, 579), chunktype=numpy.ndarray>\n",
"Coordinates: (12/21)\n",
" * time (time) datetime64[ns] 2018-03-31...\n",
" * band (band) <U12 'red' 'green' 'blue'\n",
" * x (x) float64 4.121e+05 ... 4.294e+05\n",
" * y (y) float64 4.617e+06 ... 4.605e+06\n",
" mgrs:utm_zone int64 19\n",
" proj:epsg int64 32619\n",
" ... ...\n",
" raster:bands (band) object None None None\n",
" gsd (band) object 10 10 10\n",
" common_name (band) object 'red' 'green' 'blue'\n",
" center_wavelength (band) object 0.665 0.56 0.49\n",
" full_width_half_max (band) object 0.038 0.045 0.098\n",
" epsg int64 32619