# How can I cover an AOI with a subset of image extents?

I’ve been experimenting with a problem that I have not found a friendly and simple solution for yet. Curious to know if any of you have already solved this problem.

Given a AOI and a set of image extents, how can you find a minimal number subset of images that cover? For example, the below AOI around Winnipeg must be covered by a subset of the 50 images (AOI followed by STAC query response):

The AOI:

The 50 image extents:

• I think you can use `pystac client` for this github.com/stac-utils/pystac-client Commented Oct 5, 2021 at 13:48
• I don't think `pystac-client` has this functionality built in, but it can help streamline the search and STAC object inspection. My feeling is this functionality requires a bespoke algorithm that incrementally unions Item geometries together until the unioned geometry covers the entire AOI. Commented Oct 5, 2021 at 18:38

The basic idea below is:

• find the image that covers the largest area in your AOI
• store in list
• find the next image that covers the most remaining area
• store in list
• repeat until there is no area left to cover
``````   from shapely.geometry import Polygon, MultiPolygon, box
import json
import requests

def find_largest_intersect(small_polys, large_poly):
max_area = 0
max_index = None

for i, poly in enumerate(small_polys):
intersect_poly = large_poly.intersection(poly)
if (intersect_area := intersect_poly.area) > max_area:
max_area = intersect_area
max_index = i
return max_index

# create AOI bbox
bbox = [-100.43701171875, 48.879167148960214,-95.06469726562,51.1448943093]
aoi = Polygon(box(*bbox))
area_left = aoi.area

# get STAC items
stac_api = "https://earth-search.aws.element84.com/v0/search"
limit = 50
results = requests.post(
stac_api,
data=json.dumps(
{
"bbox": bbox,
"collections": ["sentinel-s2-l2a-cogs"],
"limit": limit
}
)
)

# pull geometry from items
polygons = []
for i in results.json()["features"]:
polygons.append(Polygon(i["geometry"]["coordinates"][0]))

selected_polygons = []

counter = 0
# while there is still uncovered AOI, find the geometry that covers the most area
# and store it in `selected_polygons`
while area_left > 0 and counter <= limit:
largest_intersect_index = find_largest_intersect(polygons, aoi)
aoi = aoi.difference(polygons[largest_intersect_index])
area_left = aoi.area
selected_polygons.append(polygons.pop(largest_intersect_index))
counter += 1
print('selected', len(selected_polygons))
MultiPolygon(selected_polygons + [Polygon(box(*bbox))])
``````

Result:

• Thanks for answering @phloem! The algorithm does seem to return the correct response and I think it is fine for most uses cases with a small amount of image extents. However, finding the remaining intersecting area for all images at every iteration is an incredibly computational heavy operation. It could be improved by only recalculating the areas of the polygons that intersect the previously selected polygon. Commented Oct 7, 2021 at 12:25
• That's an excellent idea, will try Commented Oct 7, 2021 at 16:06
• I used an edge list in my POC: en.wikipedia.org/wiki/Edge_list. But there are probably other options. The complexity with doing it with the zones directly is that the intersecting zones considered should take into account only the remaining parts to be covered. I'm not sure if there's another way to avoid this constraint then by using some sort of continuous grid system on the AOI Commented Oct 7, 2021 at 16:25

I've been working on a proof of concept that finds the minimum images necessary to cover an AOI for an input set of STAC image features.

There's a lot going on and some moving pieces I have omitted out of simplicity. The general advantage of the solution is updating only the images intersecting area (or priority) based on the previous image in O(1). This is completed by creating an h3 grid of the AOI (the universe) and creating a map (key-value) where the keys are the h3 id of the universe and the values are the images containing the specific h3 id.

The core of the algorithm is laid out here:

Building the geo priority bucket:

``````package minCoverSet

import (
"fmt"
geojson "github.com/paulmach/go.geojson"
"github.com/uber/h3-go/v3"
)

type GeojsonUniverseError struct {
Value string
}

func (e *GeojsonUniverseError) Error() string {
return e.Value
}

func GetMinCoverSetGeojson(universe string, sets string) (geojson.FeatureCollection, error) {
// Get GeoPolygon of Universe
universeGeojson := []byte(universe)
ug, err  := geojson.UnmarshalFeatureCollection(universeGeojson)
if err != nil {
return geojson.FeatureCollection{}, err
}
if len(ug.Features) > 1 {
return geojson.FeatureCollection{}, &GeojsonUniverseError{Value: fmt.Sprintf("expected universe to be a single feature geojson, got: %v", ug)}
}

universePolygon, err := FeatureToGeoPolygon(*ug.Features[0])
if err != nil {
return geojson.FeatureCollection{}, err
}

h3UniverseIds := h3.Polyfill(universePolygon, 9)
values := make(map[interface{}]int, len(h3UniverseIds))
for _, id := range h3UniverseIds {
values[h3.ToString(id)] = 1
}
u := Universe{Values: values}

// Get Sets
rawSetsGeojson := []byte(sets)
setsGeojson, err  := geojson.UnmarshalFeatureCollection(rawSetsGeojson)
if err != nil {
return geojson.FeatureCollection{}, err
}

prioritySets := make([]PrioritySet, len(setsGeojson.Features))

for i, feature := range setsGeojson.Features {
featureGeoPolygon, err := FeatureToGeoPolygon(*feature)
if err != nil {
return geojson.FeatureCollection{}, err
}
prioritySets[i] = GeoSet{NameId: feature.ID, H3Ids: GetH3StringIdsFromGeoPolygon(featureGeoPolygon, 9), GeoPolygon: featureGeoPolygon, Feature: *feature}
}

// Get Min Covering sets
minSets, err := GetMinCoveringSet(prioritySets, u)

if err != nil {
return geojson.FeatureCollection{}, err
}

minSetsGeojson := geojson.NewFeatureCollection()

minFeatures := make([]geojson.Feature, len(minSets))
for i, set := range minSets {
minFeatures[i] = set.(GeoSet).Feature
}

return *minSetsGeojson, nil
}
``````

The abstract algorithm for building a minimum set:

``````package minCoverSet

import (
"container/heap"
)

type Universe struct {
Values map[interface{}]int
}

// sets are assumed to have neighbors correctly set
func GetMinCoveringSet(sets []PrioritySet, universe Universe) ([]PrioritySet, error) {
var prioritySets OrderPriorityHeap
heap.Init(&prioritySets)

universeNeighbors := make(map[interface{}][]*OrderPriority, len(sets))
for i, set := range sets {
filteredSet, err := set.Keep(universe.Values)
if err != nil {
return nil, err
}

op := OrderPriority{filteredSet, set, make(map[*OrderPriority]int, len(sets)), i} // TODO: replace nil
for v := range filteredSet.GetValues() {
if _, ok := universeNeighbors[v]; ok {
for _, n := range universeNeighbors[v] {
}
}
universeNeighbors[v] = append(universeNeighbors[v], &op) //
}
heap.Push(&prioritySets, &op)
}

// choice to iterate on workingSets and break when universe empty
// rather than while loop on universe to avoid infinite loop
// that occurs when union of workingSets != universe
var minSets []PrioritySet
size := prioritySets.Len()
for i := 0; i < size; i++ {

op := heap.Pop(&prioritySets).(OrderPriority) // Pop should update neighbors
minSets = append(minSets, op.initialSet)
op.currentSet.UpdateUniverse(&universe)
for n := range op.neighbors {
heap.Fix(&prioritySets, n.index)
}

// all values of universe found
if len(universe.Values) == 0 {
return minSets, nil
}
}

// union of sets != universe
return nil, nil

}

``````

I recorded a short demo here of a UI implementing the actual algorithm: https://youtu.be/t41Mh8Rs3k8