{ "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", "from geogif import gif, dgif" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Launch a cluster on [Coiled](https://cloud.coiled.io/) near the data." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import coiled\n", "import distributed\n", "\n", "cluster = coiled.Cluster(\n", " name=\"geogif\",\n", " backend_options={\"region\": \"us-west-1\"},\n", " n_workers=20,\n", ")\n", "client = distributed.Client(cluster)" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "bbox = (-70.05487089045347, 41.59659777453696, -69.8481905926019, 41.69921590764708)" ] }, { "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": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 1.25 s, sys: 123 ms, total: 1.37 s\n", "Wall time: 5.06 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", " )\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": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
<xarray.DataArray 'stackstac-28f5d640e4498241199302aed2dfc7aa' (time: 27,\n", " band: 3,\n", " y: 387, x: 579)>\n", "dask.array<stack, shape=(27, 3, 387, 579), dtype=float64, chunksize=(1, 2, 387, 579), chunktype=numpy.ndarray>\n", "Coordinates: (12/21)\n", " * time (time) datetime64[ns] 2017-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", " proj:epsg int64 32619\n", " mgrs:latitude_band <U1 'T'\n", " ... ...\n", " title (band) <U31 'Red (band 4) - 10m'...\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