Calculating Coverage Percentage of the AOI by an Item

This notebook demonstrates the use of pystac-client to calculate the percentage an Item’s geometry that intesects with the area of interest (AOI) specified in the search by the intersects parameter.

from pystac_client import Client

import matplotlib.pyplot as plt
from shapely.geometry import shape
from pystac_client import Client
from pystac.item import Item
from typing import Dict, Any

def intersection_percent(item: Item, aoi: Dict[str, Any]) -> float:
    '''The percentage that the Item's geometry intersects the AOI. An Item that
    completely covers the AOI has a value of 100.
    '''
    geom_item = shape(item.geometry)
    geom_aoi = shape(aoi)

    intersected_geom = geom_aoi.intersection(geom_item)

    intersection_percent = (intersected_geom.area * 100) / geom_aoi.area

    return intersection_percent


# STAC API root URL
URL = 'https://planetarycomputer.microsoft.com/api/stac/v1'

# geometry of the AOI to search over
intersects_geometry = {
    "type": "Polygon",
    "coordinates": [
        [
            [-73.21, 43.99],
            [-73.21, 47.05],
            [-70.12, 47.05],
            [-70.12, 43.99],
            [-73.21, 43.99],
        ]
    ],
}

# Create a Client and an ItemSearch representing our search
# No search operations will be performed until we call the items() method
client = Client.open(URL)
item_search = client.search(
    collections=["sentinel-2-l2a"], intersects=intersects_geometry, max_items=100
)

print([f"{intersection_percent(item, intersects_geometry):.2f}" for item in item_search.items()])

['0.32', '1.89', '1.69', '1.48', '0.35', '2.92', '6.41', '9.81', '2.96', '0.27', '0.39', '1.91', '1.69', '1.48', '0.35', '3.09', '14.93', '14.69', '14.47', '3.58', '3.59', '14.91', '14.67', '14.44', '3.15', '1.24', '7.37', '9.24', '8.94', '2.11', '0.74', '1.15', '2.60', '9.07', '5.58', '2.17', '0.06', '2.14', '9.83', '9.37', '8.81', '1.86', '0.30', '1.89', '1.69', '1.48', '0.35', '2.86', '6.35', '9.76', '2.95', '0.26', '0.39', '1.91', '1.69', '1.48', '0.35', '3.09', '14.93', '14.69', '14.47', '3.58', '3.59', '14.91', '14.67', '14.44', '3.15', '1.23', '7.34', '9.23', '8.94', '2.11', '0.77', '1.24', '2.63', '9.19', '5.69', '2.29', '0.08', '2.14', '9.83', '9.37', '8.84', '1.89', '0.29', '1.88', '1.69', '1.48', '0.35', '2.81', '6.30', '9.70', '2.94', '0.25', '0.39', '1.91', '1.69', '1.48', '0.35', '3.09']

# create a generator that filters to only those Items that intersect more than 5%
items_gt_5_percent = (
    i for i in item_search.items() if intersection_percent(i, intersects_geometry) > 5
)

# Render the AOI and Item results
# The green shape is the AOI
# The blue shapes are the Item geometries
# If there are no blue shapes, adjust the intersection percent filter above until there are

cm = plt.get_cmap('RdBu')
fig, axs = plt.subplots()
axs.set_aspect('equal', 'datalim')

for item in items_gt_5_percent:
  xs, ys = shape(item.geometry).exterior.xy
  axs.fill(xs, ys, alpha=0.5, fc='b', ec='none')

geom_intersects = shape(intersects_geometry)
xs, ys = geom_intersects.exterior.xy
axs.fill(xs, ys, alpha=0.5, fc='g', ec='none')

plt.show()