Compare the area of a coordinate polygon with the returned image from STAC

area_example.py

import matplotlib.pyplot as plt
import numpy as np
import planetary_computer as pc
import stackstac
from pystac_client import Client

import shapely.geometry
import fiona.transform

polygon = {
    "type": "Polygon",
    "coordinates": [
        [
            [-52.6762488762715, 4.08498859413025],
            [-52.6766446639803, 4.08520420586715],
            [-52.6764299428569, 4.08560163518406],
            [-52.6760341549236, 4.08538602339978],
            [-52.6762488762715, 4.08498859413025],
        ]
    ],
}


c = polygon["coordinates"][0]
c1 = [x[0] for x in c]  # first coordinate
c2 = [x[1] for x in c]
bounds = (min(c1), min(c2), max(c1), max(c2))

time_of_interest = "2015-07-01/2017-07-01"

rgb_bands = ["B04", "B03", "B02"]

resolution = 10

catalog = Client.open("https://planetarycomputer.microsoft.com/api/stac/v1")

search = catalog.search(
    collections=["sentinel-2-l2a"],
    intersects=polygon,
    datetime=time_of_interest,
    query={"eo:cloud_cover": {"lt": 20}},
)

items = search.get_all_items()
signed_items = [pc.sign(item).to_dict() for item in items]

print("epsg of items: {}".format(items[0].properties["proj:epsg"]))
geom = fiona.transform.transform_geom("epsg:4326", "epsg:32622", polygon)
coords_area = shapely.geometry.shape(geom).area

stack = stackstac.stack(
    signed_items,
    assets=rgb_bands,
    resolution=resolution,
    bounds_latlon=bounds,
    snap_bounds=False,
)

aoi = stack.compute(scheduler="single-threaded")
aoi = np.array(aoi)

img_area = aoi.shape[2] * aoi.shape[3] * resolution ** 2

print("img dims {}".format(aoi.shape))
print("Resolution per pixel in meter:{}".format(resolution))
print("Area of image: {}".format(img_area))
print("Area of polygon coordinates in m²: {}".format(coords_area))
print("Envelope area {}".format(shapely.geometry.shape(geom).envelope.area))

img = aoi[0, ...].transpose(1, 2, 0) / 2000
fig = plt.figure(figsize=(4, 4))
plt.imshow(img)
plt.show()