Access Sentinel 2 Data from AWS¶
https://registry.opendata.aws/sentinel-2-l2a-cogs/
[1]:
import dask.distributed
import folium
import folium.plugins
import geopandas as gpd
import shapely.geometry
from IPython.display import HTML, display
from pystac_client import Client
from odc.stac import configure_rio, stac_load
def convert_bounds(bbox, invert_y=False):
"""
Helper method for changing bounding box representation to leaflet notation
``(lon1, lat1, lon2, lat2) -> ((lat1, lon1), (lat2, lon2))``
"""
x1, y1, x2, y2 = bbox
if invert_y:
y1, y2 = y2, y1
return ((y1, x1), (y2, x2))
[2]:
cfg = {
"sentinel-s2-l2a-cogs": {
"assets": {
"*": {"data_type": "uint16", "nodata": 0},
"SCL": {"data_type": "uint8", "nodata": 0},
"visual": {"data_type": "uint8", "nodata": 0},
},
"aliases": {"red": "B04", "green": "B03", "blue": "B02"},
},
"*": {"warnings": "ignore"},
}
Start Dask Client¶
This step is optional, but it does improve load speed significantly. You don’t have to use Dask, as you can load data directly into memory of the notebook.
[3]:
client = dask.distributed.Client()
configure_rio(cloud_defaults=True, aws={"aws_unsigned": True}, client=client)
display(client)
Client
Client-b4582121-caaa-11ed-86c4-6045bde8e932
Connection method: Cluster object | Cluster type: distributed.LocalCluster |
Dashboard: http://127.0.0.1:8787/status |
Cluster Info
LocalCluster
0e6830c3
Dashboard: http://127.0.0.1:8787/status | Workers: 2 |
Total threads: 2 | Total memory: 6.78 GiB |
Status: running | Using processes: True |
Scheduler Info
Scheduler
Scheduler-1239b1a7-92cf-44ed-a2b0-3993ee3852c1
Comm: tcp://127.0.0.1:33975 | Workers: 2 |
Dashboard: http://127.0.0.1:8787/status | Total threads: 2 |
Started: Just now | Total memory: 6.78 GiB |
Workers
Worker: 0
Comm: tcp://127.0.0.1:37851 | Total threads: 1 |
Dashboard: http://127.0.0.1:45209/status | Memory: 3.39 GiB |
Nanny: tcp://127.0.0.1:35169 | |
Local directory: /tmp/dask-worker-space/worker-0ff_1_kd |
Worker: 1
Comm: tcp://127.0.0.1:41935 | Total threads: 1 |
Dashboard: http://127.0.0.1:41859/status | Memory: 3.39 GiB |
Nanny: tcp://127.0.0.1:36597 | |
Local directory: /tmp/dask-worker-space/worker-myzuc5jj |
Find STAC Items to Load¶
[4]:
km2deg = 1.0 / 111
x, y = (113.887, -25.843) # Center point of a query
r = 100 * km2deg
bbox = (x - r, y - r, x + r, y + r)
catalog = Client.open("https://earth-search.aws.element84.com/v0")
query = catalog.search(
collections=["sentinel-s2-l2a-cogs"], datetime="2021-09-16", limit=100, bbox=bbox
)
items = list(query.get_items())
print(f"Found: {len(items):d} datasets")
# Convert STAC items into a GeoJSON FeatureCollection
stac_json = query.get_all_items_as_dict()
Found: 9 datasets
Review Query Result¶
We’ll use GeoPandas DataFrame object to make plotting easier.
[5]:
gdf = gpd.GeoDataFrame.from_features(stac_json, "epsg:4326")
# Compute granule id from components
gdf["granule"] = (
gdf["sentinel:utm_zone"].apply(lambda x: f"{x:02d}")
+ gdf["sentinel:latitude_band"]
+ gdf["sentinel:grid_square"]
)
fig = gdf.plot(
"granule",
edgecolor="black",
categorical=True,
aspect="equal",
alpha=0.5,
figsize=(6, 12),
legend=True,
legend_kwds={"loc": "upper left", "frameon": False, "ncol": 1},
)
_ = fig.set_title("STAC Query Results")

Plot STAC Items on a Map¶
[6]:
# https://github.com/python-visualization/folium/issues/1501
from branca.element import Figure
fig = Figure(width="400px", height="500px")
map1 = folium.Map()
fig.add_child(map1)
folium.GeoJson(
shapely.geometry.box(*bbox),
style_function=lambda x: dict(fill=False, weight=1, opacity=0.7, color="olive"),
name="Query",
).add_to(map1)
gdf.explore(
"granule",
categorical=True,
tooltip=[
"granule",
"datetime",
"sentinel:data_coverage",
"eo:cloud_cover",
],
popup=True,
style_kwds=dict(fillOpacity=0.1, width=2),
name="STAC",
m=map1,
)
map1.fit_bounds(bounds=convert_bounds(gdf.unary_union.bounds))
display(fig)
Construct Dask Dataset¶
Note that even though there are 9 STAC Items on input, there is only one timeslice on output. This is because of groupby="solar_day"
. With that setting stac_load
will place all items that occured on the same day (as adjusted for the timezone) into one image plane.
[7]:
# Since we will plot it on a map we need to use `EPSG:3857` projection
crs = "epsg:3857"
zoom = 2**5 # overview level 5
xx = stac_load(
items,
bands=("red", "green", "blue"),
crs=crs,
resolution=10 * zoom,
chunks={}, # <-- use Dask
groupby="solar_day",
stac_cfg=cfg,
)
display(xx)
<xarray.Dataset> Dimensions: (y: 1099, x: 833, time: 1) Coordinates: * y (y) float64 -2.797e+06 -2.798e+06 ... -3.148e+06 -3.149e+06 * x (x) float64 1.255e+07 1.255e+07 ... 1.282e+07 1.282e+07 spatial_ref int32 3857 * time (time) datetime64[ns] 2021-09-16T02:34:44 Data variables: red (time, y, x) uint16 dask.array<chunksize=(1, 1099, 833), meta=np.ndarray> green (time, y, x) uint16 dask.array<chunksize=(1, 1099, 833), meta=np.ndarray> blue (time, y, x) uint16 dask.array<chunksize=(1, 1099, 833), meta=np.ndarray>
Note that data is not loaded yet. But we can review memory requirement. We can also check data footprint.
[8]:
xx.odc.geobox
[8]:
GeoBox
WKT
BASEGEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["unnamed",
METHOD["Popular Visualisation Pseudo Mercator",
ID["EPSG",1024]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",0,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["easting (X)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["northing (Y)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Web mapping and visualisation."],
AREA["World between 85.06°S and 85.06°N."],
BBOX[-85.06,-180,85.06,180]],
ID["EPSG",3857]]
Load data into local memory¶
[9]:
%%time
xx = xx.compute()
CPU times: user 334 ms, sys: 47.5 ms, total: 382 ms
Wall time: 5.2 s
[10]:
_ = (
xx.isel(time=0)
.to_array("band")
.plot.imshow(
col="band",
size=4,
vmin=0,
vmax=4000,
)
)

Load with bounding box¶
As you can see stac_load
returned all the data covered by STAC items returned from the query. This happens by default as stac_load
has no way of knowing what your query was. But it is possible to control what region is loaded. There are several mechanisms available, but probably simplest one is to use bbox=
parameter (compatible with stac_client
).
Let’s load a small region at native resolution to demonstrate.
[11]:
r = 6.5 * km2deg
small_bbox = (x - r, y - r, x + r, y + r)
yy = stac_load(
items,
bands=("red", "green", "blue"),
crs=crs,
resolution=10,
chunks={}, # <-- use Dask
groupby="solar_day",
stac_cfg=cfg,
bbox=small_bbox,
)
display(yy.odc.geobox)
GeoBox
WKT
BASEGEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4326]],
CONVERSION["unnamed",
METHOD["Popular Visualisation Pseudo Mercator",
ID["EPSG",1024]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",0,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["easting (X)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["northing (Y)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Web mapping and visualisation."],
AREA["World between 85.06°S and 85.06°N."],
BBOX[-85.06,-180,85.06,180]],
ID["EPSG",3857]]
[12]:
yy = yy.compute()
[13]:
_ = (
yy.isel(time=0)
.to_array("band")
.plot.imshow(
col="band",
size=4,
vmin=0,
vmax=4000,
)
)
