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:

enter image description here

The 50 image extents:


  • I think you can use pystac client for this github.com/stac-utils/pystac-client
    – GISHuman
    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.
    – phloem
    Commented Oct 5, 2021 at 18:38

2 Answers 2


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(
                "bbox": bbox,
                "collections": ["sentinel-s2-l2a-cogs"],
                "limit": limit
    # pull geometry from items
    polygons = []
    for i in results.json()["features"]:
    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
        counter += 1
    print('selected', len(selected_polygons))
    MultiPolygon(selected_polygons + [Polygon(box(*bbox))])



  • 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
    – phloem
    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 (
    geojson "github.com/paulmach/go.geojson"

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 (

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

    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)
        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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.