STAC Metadata Visualizations#
This notebook illustrates a simple way to display footprints of discovered Items after searching a STAC API, and making simple plots using Pandas and Holoviews. Only the metadata is visualized in these examples through maps and plots. The actual STAC data (i.e., Item Assets) are not accessed.
The libraries GeoPandas and hvplot are used for visualizations.
[1]:
# set pystac_client logger to DEBUG to see API calls
import logging
from pystac_client import Client
logging.basicConfig()
logger = logging.getLogger("pystac_client")
logger.setLevel(logging.DEBUG)
Define the STAC API to use, along with any custom headers (such as for authentication)
[2]:
# STAC API root URL
URL = "https://planetarycomputer.microsoft.com/api/stac/v1"
# custom headers
headers = []
cat = Client.open(URL, headers=headers)
cat
DEBUG:pystac_client.stac_api_io:GET https://planetarycomputer.microsoft.com/api/stac/v1 Headers: {'User-Agent': 'python-requests/2.34.2', 'Accept-Encoding': 'gzip, deflate', 'Accept': '*/*', 'Connection': 'keep-alive'}
[2]:
Search#
Perform a spatio-temporal search of ASTER data for a small AOI in the northern part of The Netherlands between 2000 and 2010.
[3]:
# AOI around Delfzijl, in northern Netherlands
geom = {
"type": "Polygon",
"coordinates": [
[
[6.42425537109375, 53.174765470134616],
[7.344360351562499, 53.174765470134616],
[7.344360351562499, 53.67393435835391],
[6.42425537109375, 53.67393435835391],
[6.42425537109375, 53.174765470134616],
]
],
}
# limit sets the # of items per page so we can see multiple pages getting fetched
search = cat.search(
max_items=50,
collections="aster-l1t",
intersects=geom,
datetime="2000-01-01/2010-12-31",
)
# retrieve the items as dictionaries, rather than Item objects
items = list(search.items_as_dicts())
len(items)
DEBUG:pystac_client.stac_api_io:GET https://planetarycomputer.microsoft.com/api/stac/v1/ Headers: {'User-Agent': 'python-requests/2.34.2', 'Accept-Encoding': 'gzip, deflate', 'Accept': '*/*', 'Connection': 'keep-alive'}
DEBUG:pystac_client.stac_api_io:POST https://planetarycomputer.microsoft.com/api/stac/v1/search Headers: {'User-Agent': 'python-requests/2.34.2', 'Accept-Encoding': 'gzip, deflate', 'Accept': '*/*', 'Connection': 'keep-alive', 'Content-Length': '341', 'Content-Type': 'application/json'} Payload: {"datetime": "2000-01-01T00:00:00Z/2010-12-31T23:59:59Z", "collections": ["aster-l1t"], "intersects": {"type": "Polygon", "coordinates": [[[6.42425537109375, 53.174765470134616], [7.344360351562499, 53.174765470134616], [7.344360351562499, 53.67393435835391], [6.42425537109375, 53.67393435835391], [6.42425537109375, 53.174765470134616]]]}}
[3]:
50
GeoPandas#
A GeoDataFrame is constructed from the AOI geometry, in order to make visualizations easy.
The STAC Items, which are a GeoJSON FeatureCollection can be converted to a GeoDataFrame.
[4]:
from copy import deepcopy
import geopandas as gpd
import pandas as pd
from shapely.geometry import shape
# convert a list of STAC Items into a GeoDataFrame
def items_to_geodataframe(items):
_items = []
for i in items:
_i = deepcopy(i)
_i["geometry"] = shape(_i["geometry"])
_items.append(_i)
gdf = gpd.GeoDataFrame(pd.json_normalize(_items))
for field in ["properties.datetime", "properties.created", "properties.updated"]:
if field in gdf:
gdf[field] = pd.to_datetime(gdf[field])
gdf.set_index("properties.datetime", inplace=True)
return gdf
[5]:
# convert geometry to a GeoDataFrame
aoi_gdf = gpd.GeoDataFrame([{"geometry": shape(geom)}])
aoi_gdf
[5]:
| geometry | |
|---|---|
| 0 | POLYGON ((6.42426 53.17477, 7.34436 53.17477, ... |
[6]:
# convert found items to a GeoDataFrame
items_gdf = items_to_geodataframe(items)
items_gdf.head()
[6]:
| id | bbox | type | links | geometry | collection | stac_extensions | stac_version | assets.TIR.href | assets.TIR.type | ... | assets.vnir-browse.href | assets.vnir-browse.type | assets.vnir-browse.roles | assets.vnir-browse.title | assets.vnir-browse.description | assets.qa-txt.href | assets.qa-txt.type | assets.qa-txt.roles | assets.qa-txt.title | assets.qa-txt.description | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| properties.datetime | |||||||||||||||||||||
| 2006-12-14 21:25:39.910000+00:00 | AST_L1T_00312142006212539_20150517105406 | [5.8893619, 53.006463, 7.079664, 53.6760341] | Feature | [{'rel': 'collection', 'type': 'application/js... | POLYGON ((6.83853 53.67589, 6.84669 53.67603, ... | aster-l1t | [https://stac-extensions.github.io/eo/v1.0.0/s... | 1.0.0 | https://astersa.blob.core.windows.net/aster/im... | image/tiff; application=geotiff; profile=cloud... | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 2006-09-19 10:50:24.192000+00:00 | AST_L1T_00309192006105024_20150516061550 | [6.2313306, 52.613738, 7.4219194, 53.2988665] | Feature | [{'rel': 'collection', 'type': 'application/js... | POLYGON ((7.42138 53.15465, 7.42093 53.15377, ... | aster-l1t | [https://stac-extensions.github.io/eo/v1.0.0/s... | 1.0.0 | https://astersa.blob.core.windows.net/aster/im... | image/tiff; application=geotiff; profile=cloud... | ... | https://astersa.blob.core.windows.net/aster/im... | image/jpeg | [thumbnail] | VNIR browse file | Standalone reduced resolution VNIR | NaN | NaN | NaN | NaN | NaN |
| 2006-09-19 10:50:15.355000+00:00 | AST_L1T_00309192006105015_20150516061540 | [6.4862887, 53.1318137, 7.6929281, 53.8181412] | Feature | [{'rel': 'collection', 'type': 'application/js... | POLYGON ((7.69274 53.67237, 7.69231 53.67153, ... | aster-l1t | [https://stac-extensions.github.io/eo/v1.0.0/s... | 1.0.0 | https://astersa.blob.core.windows.net/aster/im... | image/tiff; application=geotiff; profile=cloud... | ... | https://astersa.blob.core.windows.net/aster/im... | image/jpeg | [thumbnail] | VNIR browse file | Standalone reduced resolution VNIR | https://astersa.blob.core.windows.net/aster/im... | text/plain | [metadata] | QA browse file | Geometric quality assessment report. |
| 2006-07-17 10:50:29.170000+00:00 | AST_L1T_00307172006105029_20150515082245 | [5.4068552, 53.3010121, 6.5997335, 53.9797445] | Feature | [{'rel': 'collection', 'type': 'application/js... | POLYGON ((6.095 53.91815, 6.11773 53.91475, 6.... | aster-l1t | [https://stac-extensions.github.io/eo/v1.0.0/s... | 1.0.0 | https://astersa.blob.core.windows.net/aster/im... | image/tiff; application=geotiff; profile=cloud... | ... | https://astersa.blob.core.windows.net/aster/im... | image/jpeg | [thumbnail] | VNIR browse file | Standalone reduced resolution VNIR | https://astersa.blob.core.windows.net/aster/im... | text/plain | [metadata] | QA browse file | Geometric quality assessment report. |
| 2006-07-03 10:38:13.100000+00:00 | AST_L1T_00307032006103813_20150515031748 | [7.2238549, 52.9288706, 8.4058086, 53.6021029] | Feature | [{'rel': 'collection', 'type': 'application/js... | POLYGON ((8.37763 53.47724, 8.40581 53.47338, ... | aster-l1t | [https://stac-extensions.github.io/eo/v1.0.0/s... | 1.0.0 | https://astersa.blob.core.windows.net/aster/im... | image/tiff; application=geotiff; profile=cloud... | ... | https://astersa.blob.core.windows.net/aster/im... | image/jpeg | [thumbnail] | VNIR browse file | Standalone reduced resolution VNIR | https://astersa.blob.core.windows.net/aster/im... | text/plain | [metadata] | QA browse file | Geometric quality assessment report. |
5 rows × 76 columns
Plot Geometries on a Map#
Holoviews is used to display geometries on a map by using hvplot. The The * Holoviews operator to overlay two plots
[7]:
import hvplot.pandas # noqa: F401
# plot polygons on a map with background tiles.
def plot_polygons(data, *args, **kwargs):
return data.hvplot.polygons(
*args,
geo=True,
projection="GOOGLE_MERCATOR",
xaxis=None,
yaxis=None,
frame_width=600,
frame_height=600,
fill_alpha=0,
line_width=4,
**kwargs,
)
[8]:
plot_polygons(aoi_gdf, tiles="OSM", line_color="red")
[8]:
[9]:
plot_polygons(items_gdf, tiles="OSM")
[9]:
[10]:
plot_polygons(items_gdf, tiles="OSM") * plot_polygons(aoi_gdf, line_color="red")
[10]:
Line Plots#
Numeric STAC metadata can also be plotted, most often this will be plotted vs the Item datetime.
[11]:
items_df = pd.DataFrame(items_gdf)
plot_fields = [
"properties.eo:cloud_cover",
"properties.view:sun_azimuth",
"properties.view:sun_elevation",
]
items_df[plot_fields].hvplot(height=500, width=800).opts(legend_position="top_right")
[11]: