Examples¶

We’ll watch coastline movement from longshore drift on Cape Cod, Massachusetts, USA.

We’ll make quarterly, cloud-free-ish composites of Sentinel-2 imagery from AWS. Then turn them into GIFs.

[1]:

import pystac_client
import stackstac
import dask.array as da
from geogif import gif, dgif

Launch a cluster on Coiled near the data.

[2]:

import coiled
import distributed

cluster = coiled.Cluster(
    name="geogif",
    backend_options={"region": "us-west-1"},
    n_workers=20,
)
client = distributed.Client(cluster)

[2]:

bbox = (-70.05487089045347, 41.59659777453696, -69.8481905926019, 41.69921590764708)

Search for STAC items overlapping our area of interest, from AWS’s open data catalog.

[3]:

%%time
items = (
    pystac_client.Client.open("https://earth-search.aws.element84.com/v1")
    .search(
        bbox=bbox,
        collections=["sentinel-2-l2a"]
    )
).item_collection()

CPU times: user 1.25 s, sys: 123 ms, total: 1.37 s
Wall time: 5.06 s

Make quarterly median composites¶

Turn STAC items into xarray as a temporal stack, using stackstac.

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.

Each 3-month composite will be one frame in our animation. We end up with 18 frames.

[4]:

stack = stackstac.stack(items, bounds_latlon=bbox, resolution=30)

scl = stack.sel(band=["scl"])
# Sentinel-2 Scene Classification Map: nodata, saturated/defective, dark, cloud shadow, cloud med. prob., cloud high prob., cirrus
invalid = da.isin(scl, [0, 1, 2, 3, 8, 9, 10])
valid = stack.where(~invalid)

rgb = valid.sel(band=["red", "green", "blue"])

quarterly = rgb.resample(time="Q").median()
quarterly

[4]:

<xarray.DataArray 'stackstac-28f5d640e4498241199302aed2dfc7aa' (time: 27,
                                                                band: 3,
                                                                y: 387, x: 579)>
dask.array<stack, shape=(27, 3, 387, 579), dtype=float64, chunksize=(1, 2, 387, 579), chunktype=numpy.ndarray>
Coordinates: (12/21)
  * time                                     (time) datetime64[ns] 2017-03-31...
  * band                                     (band) <U12 'red' 'green' 'blue'
  * x                                        (x) float64 4.121e+05 ... 4.294e+05
  * y                                        (y) float64 4.617e+06 ... 4.605e+06
    proj:epsg                                int64 32619
    mgrs:latitude_band                       <U1 'T'
    ...                                       ...
    title                                    (band) <U31 'Red (band 4) - 10m'...
    gsd                                      (band) object 10 10 10
    common_name                              (band) object 'red' 'green' 'blue'
    center_wavelength                        (band) object 0.665 0.56 0.49
    full_width_half_max                      (band) object 0.038 0.045 0.098
    epsg                                     int64 32619

Use .persist() to start the computation (only takes ~45sec), then hold the pre-computed composites in the cluster’s memory. We’re going to reuse the same composites a few times, so it makes sense to persist it.

[18]:

ts = quarterly.persist()

Quick filmstrip of all the frames, to see what we’re dealing with¶

This is actually much slower than computing a GIF itself, since we have to download the full-bit-depth data from the cluster. The data transfer usually takes longer than the computation, depending on your internet speeds.

We’re just showing it here for comparison to the animations.

Frame #2 looks so bad, let’s just plan on skipping it…

[8]:

ts_local = ts.compute()
ts_local.plot.imshow(col="time", rgb="band", col_wrap=5, robust=True)

[8]:

<xarray.plot.facetgrid.FacetGrid at 0x13e97b790>

Clean up NaNs¶

Forward-fill, the backwards-fill the NaNs, to make a smoother-looking animation. (Note that this requires Bottleneck.)

Also drop the first 2 frames.

[19]:

cleaned = ts[2:].ffill("time").bfill("time")

Make GIFs!¶

Defaults¶

[24]:

dgif(cleaned).compute()

[24]:

No timestamp¶

[25]:

dgif(cleaned, date_format=None).compute()

[25]:

S L O W M O T I O N¶

[26]:

dgif(cleaned, fps=10).compute()

[26]:

`gif` (non-Dask)¶

gif and dgif behave the same. dgif(x).compute() is usually faster than gif(x.compute()).

[17]:

gif(ts[2:])

[17]:

`gif` to file¶

[19]:

gif(ts[2:], to="capecod.gif")

`dgif` to bytes¶

[20]:

gif_bytes = dgif(cleaned, bytes=True).compute()

[21]:

with open("capecod2.gif", "wb") as f:
    f.write(gif_bytes)

Colormapped (single-band) GIF¶

We’ll demonstrate single-band colormapped GIFs with a Normalized Difference Vegetation Index (NDVI).

(There isn’t much vegetation to look at here, so this isn’t very interesting.)

[28]:

nir, red = valid.sel(band="red"), valid.sel(band="nir")
ndvi = (nir - red) / (nir + red)
quarterly_ndvi = ndvi.resample(time="Q").median()
ndvi_ts = quarterly_ndvi[2:].compute()

Defaults to viridis

[29]:

gif(ndvi_ts)

[29]:

Pick a colormap¶

[30]:

gif(ndvi_ts, cmap="YlGn")

[30]:

Custom date format and position¶

[31]:

gif(ndvi_ts, cmap="YlGn", date_format="NDVI: %Y", date_position="lr")

[31]:

Custom date color and background¶

[32]:

gif(ndvi_ts, cmap="YlGn", date_color=(0, 10, 0), date_bg=None, date_position="lr")

[32]: