Find and retrieve Sentinel-2 imagery¶

This notebook demonstrates how to:

intersect Sentinel-2 scenes with a flight track using BigQuery
retrieve and visualized Sentinel-2 scenes

This notebook relies on optional dependencies in the sat group (pip install pycontrails[sat]).

[1]:

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pyproj
from matplotlib import dates

from pycontrails.core import Flight
from pycontrails.datalib import sentinel

Intersect Sentinel-2 scenes with IAGOS flight¶

Intersecting Sentinel-2 scenes with flight tracks requires a Google Cloud account and project with access to the BigQuery API. Users without BigQuery API access can skip this section.

[2]:

df = pd.read_csv("data/iagos-flight-sentinel.csv", parse_dates=["time"])
# IAGOS in-situ data has a small time delay. For most applications, this small
# error can be ignored. For comparison with LEO satellite data, however, it is
# important to correct for this delay. The delay for this flight is approximately
# 9 seconds (this can be determined by comparing with ADS-B).
df["time"] += pd.Timedelta(seconds=9)
flight = Flight(df)

[3]:

scenes = sentinel.intersect(flight)
scenes

[3]:

	base_url	granule_id	sensing_time	source_url
5	gs://gcp-public-data-sentinel-2/tiles/30/U/WG/...	L1C_T30UWG_A014649_20191226T112436	2019-12-26 11:25:47.117000+00:00	https://storage.googleapis.com/gcp-public-data...

Retrieve and visualize Sentinel-2 scene¶

This section retrieves Sentinel-2 imagery from a Google Cloud public dataset. Retrieval uses anonymous access, so this section of the notebook can be run by users without a Google Cloud account.

The base_url, granule_id, and sensing_time variables are hard-coded so this section can be run independently of the previous section.

[4]:

base_url = "gs://gcp-public-data-sentinel-2/tiles/30/U/WG/S2B_MSIL1C_20191226T112359_N0208_R037_T30UWG_20191226T115227.SAFE"
granule_id = "L1C_T30UWG_A014649_20191226T112436"
sensing_time = pd.Timestamp("2019-12-26 11:25:47.117000+00:00")

The IAGOS flight is in an ISSR at the time of imagery acquisition:

[5]:

plt.plot(df["time"], df["rhi"], "k-")
plt.ylabel("RHi (nondim)")
plt.gca().xaxis.set_major_formatter(dates.DateFormatter("%m/%d\n%H:%M"))
plt.gca().axhline(y=1, color="gray", ls="--")
plt.gca().axvline(x=sensing_time, color="gray", ls="--");

Naive Overlay¶

We can project the 1-dimensional flight track into the coordinate reference system of the Sentinel-2 imagery and overlay it on the true-color image. The BigQuery table used to intersect flight tracks with Sentinel-2 images provides a single nominal sensing time for each image, and we can approximately locate the aircraft by identifying the point along the projected flight track that corresponds to the nominal sensing time. This point is shown below.

[8]:

# The crs artifact can be used to project flight latitude and longitude into image coordinates
transformer = pyproj.Transformer.from_crs("EPSG:4326", crs, always_xy=True)

flight_x, flight_y = transformer.transform(df["longitude"], df["latitude"])
flight_z = flight.altitude  # needed below

# Convert time to ns since epoch for interpolation
x0 = np.interp(sensing_time.value, flight["time"].astype(int), flight_x)
y0 = np.interp(sensing_time.value, flight["time"].astype(int), flight_y)

[9]:

# seams in the image are from different detectors in the Sentinel-2 sensor
plt.figure(figsize=(10, 10))
plt.imshow(rgb, extent=extent)
plt.plot(x0, y0, "ro")
plt.xlabel("Easting (m)")
plt.ylabel("Northing (m)");

If we zoom in, we can see the aircraft and contrail! We also see that the expected location of the aircraft is somewhat different from what the satellite observed.

[10]:

plt.figure(figsize=(10, 10))
plt.imshow(rgb, extent=extent)
plt.plot(x0, y0, "ro")
plt.xlim([x0 - 4e3, x0 + 4e3])
plt.ylim([y0 - 4e3, y0 + 4e3])
plt.xlabel("Easting (m)")
plt.ylabel("Northing (m)");

Corrected Overlay¶

The naive overlay has two main sources of error:

The Sentinel-2 imagery is geolocated at the surface, whereas the aircraft is at an altitude of roughly 10 km. This introduces a parallax error when the aircraft is not directly above the ground point being imaged.
The sensing time varies across the image. Instead of interpolating the flight track to the sensing time at the image center, we must interpolate the flight track to the sensing time at each pixel.

The Sentinel-2 metadata provides the information needed to correct for these two errors. The corrected overlay is shown below.

[11]:

x1, y1, t1 = handler.colocate_flight(flight)

As shown below, applying these corrections improves the agreement between the expected and observed aircraft position. Some sources of error still remain:

GPS position uncertainty: ADS-B reports aircraft positions with typical errors on the order of tens of meters.
Timestamp uncertainty: the ADS-B specification does not include a timestamp field. Receivers assign timestamps, introducing two errors: (1) receiver time may differ from the aircraft’s GPS time, and (2) timestamps usually have only 1-second resolution, which is coarse relative to satellite scan times.
Surface-based geolocation: Sentinel-2 imagery is geolocated to the Earth’s surface, assumed to be sea level in the current pycontrails correction code. This introduces parallax error that grows with surface elevation.

[12]:

plt.figure(figsize=(10, 10))
plt.imshow(rgb, extent=extent)
plt.plot(x1, y1, "ro")
plt.xlim([x1 - 4e3, x1 + 4e3])
plt.ylim([y1 - 4e3, y1 + 4e3])
plt.xlabel("Easting (m)")
plt.ylabel("Northing (m)");

Cirrus band¶

Sentinel-2 does not have true thermal bands (i.e., bands where top-of-atmosphere radiances are dominated by emission rather than reflection). However, aircraft and contrails often show up clearly in the cirrus band (B10), which is designed to isolate high clouds.

This section also demonstrates plotting with the xarray.Dataset returned by Sentinel.get.

[13]:

handler = sentinel.Sentinel(base_url, granule_id, bands=["B10"])
ds = handler.get()

[14]:

plt.figure(figsize=(12, 10))
ds["B10"].plot(vmin=0.3, vmax=0.5, cmap="gray")
plt.gca().set_aspect("equal", "box")
plt.plot(x1, y1, "ro")
plt.xlim([x1 - 4e3, x1 + 4e3])
plt.ylim([y1 - 4e3, y1 + 4e3])
plt.xlabel("Easting (m)")
plt.ylabel("Northing (m)");