Data discovery

The first step in many workflows involves discovering individual spatial data files covering the space, time, and variables of interest. Here we use a STAC Catalog API to recover a list of candidate data. We dig deeper into how this works and what it returns in later recipes. This example searches for images in a lon-lat bounding box from a collection of Cloud-Optimized-GeoTIFF (COG) images taken by Sentinel2 satellite mission. This function will not download any imagery, it merely gives us a list of metadata about available images, including the access URLs.

box = [-122.51, 37.71, -122.36, 37.81]
items = (
  Client.
  open("https://earth-search.aws.element84.com/v1").
  search(
    collections = ['sentinel-2-l2a'],
    bbox = box,
    datetime = "2022-06-01/2022-08-01",
    query={"eo:cloud_cover": {"lt": 20}}).
  item_collection()
)

We pass this list of images to a high-level utilty (gdalcubes in R, odc.stac in python) that will do all of the heavy lifting. Using the URLs and metadata provided by STAC, these functions can extract only our data of interest (given by the bounding box) without downloading unnecessary regions or bands. While streaming the data, these functions will also reproject it into the desired coordinate reference system – (an often costly operation to perform in R) and can potentially resample or aggregate the data to a desired spatial resolution. (The R code will also resample from images in overlapping areas to replace pixels masked by clouds)

data = odc.stac.load(
    items,
    crs="EPSG:32610",
    bands=["nir08", "red"],
    bbox=box
)
# For some reason, requesting reprojection in EPSG:4326 gives only empty values

We can do arbitrary calculations on this data as well. Here we calculate NDVI, a widely used measure of greenness that can be used to determine vegetation cover. (Note that the R example uses lazy evaluation, and can thus perform these calculations while streaming)

ndvi = (((data.nir08 - data.red) / (data.red + data.nir08)).
        resample(time="MS").
        median("time", keep_attrs=True).
        compute()
)

# mask out bad pixels
ndvi = ndvi.where(ndvi <= 1)

And we plot the result. The long rectangle of Golden Gate Park is clearly visible in the North-West.

import matplotlib as plt
cmap = plt.colormaps.get_cmap('viridis')  # viridis is the default colormap for imshow
cmap.set_bad(color='black')

ndvi.plot.imshow(row="time", cmap=cmap, add_colorbar=False, size=4)

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

Zonal statistics

In addition to large scale raster data such as satellite imagery, the analysis of vector shapes such as polygons showing administrative regions is a central component of spatial analysis, and particularly important to spatial social sciences. The red-lined areas of the 1930s are one example of spatial vectors. One common operation is to summarise the values of all pixels falling within a given polygon, e.g. computing the average greenness (NDVI)

# make sure we are in the same projection.
# zonal_stats method understands paths but not xarray
(   ndvi.
    rio.reproject("EPSG:4326").
    rio.to_raster(raster_path="ndvi.tif", driver="COG")   
)

sf_url = "/vsicurl/https://dsl.richmond.edu/panorama/redlining/static/citiesData/CASanFrancisco1937/geojson.json"
mean_ndvi = zonal_stats(sf_url, "ndvi.tif", stats="mean")

We plot the underlying NDVI as well as the average NDVI of each polygon, along with it’s textual grade. Note that “A” grades tend to be darkest green (high NDVI) while “D” grades are frequently the least green. (Regions not zoned for housing at the time of the 1937 housing assessment are not displayed as polygons.)

sf = gpd.read_file(sf_url)
sf["ndvi"] = [x["mean"] for x in mean_ndvi]
sf.plot(column="ndvi", legend=True)

<Axes: >

Are historically redlined areas still less green?

import geopolars as gpl
import polars as pl

(gpl.
  from_geopandas(sf).
  group_by("grade").
  agg(pl.col("ndvi").mean()).
  sort("grade")
)