AWS COG And STAC example

aws_cog_and_stac_example.py

from pyproj import Transformer
import rasterio
from rasterio.windows import from_bounds

def get_window(src, lat_north, lon_west, lat_south, lon_east):
    """ transfrom crs and get window """
    Transf = Transformer.from_crs("epsg:4326", src.crs) 
    lat_north, lon_west = Transf.transform(lat_north, lon_west)
    lat_south, lon_east = Transf.transform(lat_south, lon_east)         

    x_top, y_top = src.index( lat_north, lon_west )
    x_bottom, y_bottom = src.index( lat_south, lon_east )
    
    # Define window in RasterIO
    window = Window.from_slices( ( x_top, x_bottom ), ( y_top, y_bottom ), boundless=True ) 
    return window

factor = 1e-4
# https://www.matecdev.com/posts/landsat-sentinel-aws-s3-python.html
def getSubset(red_path, green_path, blue_path, filename, bbox):
    with rasterio.open(red_path) as red_src, \
        rasterio.open(green_path) as green_src, \
        rasterio.open(blue_path) as blue_src:

        # Download data specified by window
        srcs = [red_src, green_src, blue_src]
        lon_west, lat_south, lon_east, lat_north = bbox
        window = get_window(red_src, lat_north, lon_west, lat_south, lon_east)
        # Actual HTTP range request
        bands = [src.read(1, window=window) * factor for src in srcs]

        transform = rasterio.windows.transform(window, red_src.transform)
        # Save
        kwargs = red_src.meta.copy()
        kwargs.update({
            'height': bands[0].shape[0],
            'width': bands[0].shape[1],
            'count': len(srcs),
            'dtype': bands[0].dtype,
            'transform': transform
        })

        with rasterio.open(f"./{filename}.tif", "w", **kwargs) as new_dataset:
            for i, band in enumerate(bands):
                new_dataset.write(band, i+1)
                
                
# Define the geographical area of interest
center_latitude = 35.308431
center_longitude = 139.319674
deg_range = 0.01 # degree
cloud_cover_thresh = 20 # percentage
datetime = "2023-05-01T00:00:00Z/2023-05-31T23:59:59Z"

bbox = [
    center_longitude - deg_range, 
    center_latitude - deg_range, 
    center_longitude + deg_range, 
    center_latitude + deg_range
]
search = Search(
    bbox=bbox,
    datetime=datetime,
    collections=["sentinel-2-l2a"],
    url="https://earth-search.aws.element84.com/v1",
    limit=100,
)

items = search.items()

for i, item in enumerate(items):
    cloud_cover = item["eo:cloud_cover"]
    if cloud_cover < cloud_cover_thresh:
        print(i, item["datetime"], cloud_cover)
        
# Use low cloud coverage.
index = 1
item = items[index]


red_path = item.assets["red"]["href"]
green_path = item.assets["green"]["href"]
blue_path = item.assets["blue"]["href"]

with rasterio.open(href) as src:
    getSubset(red_path, green_path, blue_path, str(item), bbox)