Accessing MODIS temperature data with the Planetary Computer STAC API¶

The planetary computer hosts three temperature-related MODIS 6.1 products:

Land Surface Temperature/Emissivity Daily (11A1)
Land Surface Temperature/Emissivity 8-Day (11A2)
Land Surface Temperature/3-Band Emissivity 8-Day (21A2)

For more information about the products themselves, check out the User Guides at the bottom of this document.

Environment setup¶

This notebook works with or without an API key, but you will be given more permissive access to the data with an API key. The Planetary Computer Hub is pre-configured to use your API key.

In [1]:

import odc.stac
import planetary_computer
import pystac_client
import rich.table

Data access¶

The datasets hosted by the Planetary Computer are available from Azure Blob Storage. We'll use pystac-client to search the Planetary Computer's STAC API for the subset of the data that we care about, and then we'll load the data directly from Azure Blob Storage. We'll specify a modifier so that we can access the data stored in the Planetary Computer's private Blob Storage Containers. See Reading from the STAC API and Using tokens for data access for more.

In [2]:

catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1",
    modifier=planetary_computer.sign_inplace,
)

Query for available data¶

MODIS is a global dataset with a variety of products available within each larger category (vegetation, snow, fire, temperature, and reflectance). The MODIS group contains a complete listing of available collections. Each collection's format followsmodis-{product}-061, where product is the MODIS product id. The -061 suffix indicates that all of the MODIS collections are part of the MODIS 6.1 update.

Let's access Land Surface Temperature/Emissivity Daily (11A1) data over Boise, Idaho in 2021. We'll get four images for the midseasonal months: March, June, September, and December.

In [3]:

# Boise, Idaho
latitude = 43.6
longitude = -116.2
buffer = 1
bbox = [longitude - buffer, latitude - buffer, longitude + buffer, latitude + buffer]
year = "2021"
months = {
    "March": "03",
    "June": "06",
    "September": "09",
    "December": "12",
}
items = dict()

# Fetch the collection of interest and print available items
for name, number in months.items():
    datetime = f"{year}-{number}"
    search = catalog.search(
        collections=["modis-11A1-061"],
        bbox=bbox,
        datetime=datetime,
    )
    items[name] = search.get_all_items()[0]

print(items)

{'March': <Item id=MYD11A1.A2021090.h09v04.061.2021092201848>, 'June': <Item id=MYD11A1.A2021181.h09v04.061.2021183191149>, 'September': <Item id=MYD11A1.A2021273.h09v04.061.2021277222441>, 'December': <Item id=MYD11A1.A2021365.h09v04.061.2022004014918>}

Available assets¶

Each item has several available assets, including the original HDF file and a Cloud-optimized GeoTIFF of each subdataset.

In [4]:

t = rich.table.Table("Key", "Title")
for key, asset in items["March"].assets.items():
    t.add_row(key, asset.title)
t

Out[4]:

┏━━━━━━━━━━━━━━━━━━┳━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┓
┃ Key              ┃ Title                                                               ┃
┡━━━━━━━━━━━━━━━━━━╇━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┩
│ hdf              │ Source data containing all bands                                    │
│ QC_Day           │ Quality control for daytime LST and emissivity                      │
│ Emis_31          │ Band 31 emissivity                                                  │
│ Emis_32          │ Band 32 emissivity                                                  │
│ QC_Night         │ Quality control for nighttime LST and emissivity                    │
│ metadata         │ Federal Geographic Data Committee (FGDC) Metadata                   │
│ LST_Day_1km      │ Daily daytime 1km grid Land-surface Temperature                     │
│ Clear_day_cov    │ Day clear-sky coverage                                              │
│ Day_view_angl    │ View zenith angle of daytime Landsurface Temperature                │
│ Day_view_time    │ (local solar) Time of daytime Land-surface Temperature observation  │
│ LST_Night_1km    │ Daily nighttime 1km grid Land-surface Temperature                   │
│ Clear_night_cov  │ Night clear-sky coverage                                            │
│ Night_view_angl  │ View zenith angle of nighttime Landsurface Temperature              │
│ Night_view_time  │ (local solar) Time of nighttime Landsurface Temperature observation │
│ tilejson         │ TileJSON with default rendering                                     │
│ rendered_preview │ Rendered preview                                                    │
└──────────────────┴─────────────────────────────────────────────────────────────────────┘

Loading the data¶

For this example, we'll visualize the temperature data over Boise, Idaho. Let's grab each fire mask cover COG and load them into an xarray using odc-stac. The MODIS coordinate reference system is a sinusoidal grid, which means that views in a naïve XY raster look skewed. For visualization purposes, we reproject to a spherical Mercator projection for intuitive, north-up visualization.

In [5]:

data = odc.stac.load(
    items.values(),
    crs="EPSG:3857",
    bands="LST_Day_1km",
    resolution=500,
    bbox=bbox,
)

raster = items["March"].assets["LST_Day_1km"].extra_fields["raster:bands"]
data = data["LST_Day_1km"] * raster[0]["scale"]
data

Displaying the data¶

Let's display the temperature for each month.

In [6]:

g = data.plot.imshow(cmap="magma", col="time", vmin=250, vmax=325, size=4)
datetimes = data.time.to_pandas().dt.strftime("%B")

for ax, datetime in zip(g.axes.flat, datetimes):
    ax.set_title(datetime)

Change detection¶

Now let's see how the temperature changes between seasons.

In [7]:

change = data.diff(dim="time").assign_coords(time=["Winter", "Spring", "Summer"])
change.plot.imshow(cmap="coolwarm_r", col="time");