Minimal example demonstrating how to pull Sentinel-2 data for a specific location.

sentinel2_example.py

# Heavily inspired by: https://geospatial.101workbook.org/ImportingData/ImportingImages.html.

import geopandas as gpd
import numpy as np
import pandas as pd
import requests
import stackstac

from PIL import Image
from satsearch import Search

df = pd.DataFrame({"latitude": [32.60289673471132], "longitude": [-85.48890553734162]})
gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df["longitude"], df["latitude"]))
gdf.set_crs("EPSG:4326", inplace=True)
wkt_url = "https://spatialreference.org/ref/esri/usa-contiguous-albers-equal-area-conic/prettywkt/"
wkt = requests.get(wkt_url).text
gdf.to_crs(wkt, inplace=True)
radius_meters = 400
# Defines a circle with a 400 m radius centered at the GPS coordinates.
gdf["geometry"] = gdf["geometry"].buffer(radius_meters)

# Query the catalog using the SAT-API: http://sat-utils.github.io/sat-api/.
sent2_proj = 4326
search_args = {
    # Endpoint for accessing Sentinel-2 data hosted on AWS through their Open Data
    # initiative.
    "url": "https://earth-search.aws.element84.com/v0",
    "collections": ["sentinel-s2-l2a-cogs"],
    "datetime": "2021-09-01/2021-09-15",
    "bbox": list(gdf.to_crs(sent2_proj).total_bounds),
}
# STAC Items are described here: https://github.com/radiantearth/stac-spec/blob/master/item-spec/item-spec.md.
stac_items = Search.search(**search_args).items()

# Transform the assets into an xarray object, reading all the bands. A description of
# the various bands can be found at: https://gisgeography.com/sentinel-2-bands-combinations/.
bands = [f"B{str(idx).zfill(2)}" for idx in range(1, 13)]
da = stackstac.stack(stac_items.geojson()["features"], assets=bands)

# Get the bounds of the region of interest (ROI) using the projection found in the
# DataArray. This bounding box is used to clip the Sentinel-2 Level-2A tiles to the
# extent of the ROI.
(minx, miny, maxx, maxy) = gdf.to_crs(int(da["epsg"])).total_bounds
da_sub = da.sel(x=slice(minx, maxx), y=slice(maxy, miny))

# Convert the DataArray to a Dataset.
data_name = "data"
da_sub = da_sub.to_dataset(name=data_name)
da_sub = da_sub[data_name]

# Download the data from AWS. Takes a few seconds.
da_sub = da_sub.compute()

# See: https://www.mdpi.com/2072-4292/9/6/584/htm.
# "The surface reflectance values are coded in JPEG2000 with the same quantification
# value of 10,000 as for Level-1C products, i.e., a factor of 1/10,000 needs to be
# applied to Level-2A digital numbers (DN) to retrieve physical surface reflectance
# values."
max_val = 10000

# Loop through the dates with imagery. Notice that 2021-09-01 and 2021-09-11 are cloudy
# days.
for dt in da_sub["time"].values:
    print(dt)

    # Get the data for a single date and only the RGB bands.
    da_vis = da_sub.sel(time=str(dt), band=["B04", "B03", "B02"])

    # Convert the data into an image.
    img_arr = da_vis.squeeze().values.copy()
    img_arr /= max_val
    img_arr = img_arr.clip(0.0, 1.0)
    img_arr = np.transpose(img_arr, (1, 2, 0))
    img = Image.fromarray(np.array(255 * img_arr, dtype="uint8"))
    img.resize((512, 512), resample=Image.LANCZOS).show()

    # Convert the data for each of the other bands into an image.
    for band in range(5, 13):
        try:
            da_band = da_sub.sel(time=str(dt), band=f"B{str(band).zfill(2)}")
        except KeyError:
            continue

        img_arr = da_band.squeeze().values.copy()
        img_arr /= max_val
        img_arr = img_arr.clip(0.0, 1.0)
        img = Image.fromarray(np.array(255 * img_arr, dtype="uint8"))
        img.resize((512, 512), resample=Image.LANCZOS).show()