gjoseph92 / stackstac

Turn a STAC catalog into a dask-based xarray
https://stackstac.readthedocs.io
MIT License
248 stars 49 forks source link

Adapting Basic Example for L1C Sentinel-2 Products #211

Closed J6767 closed 1 year ago

J6767 commented 1 year ago

The basic example runs great for me. Now, I am trying to adapt the basic example here but for Sentinel-2 L1C products (which are needed for some specific atmospheric processing). I cannot use L2 for this.

from pystac_client import Client
import stackstac
import pyproj
import dask.diagnostics
from rasterio.session import AWSSession
import rasterio as rio

client = Client.open("https://earth-search.aws.element84.com/v1")
lon, lat = -105.78, 35.79

search = client.search(
    max_items=10,
    collections=['sentinel-2-l1c'],
    # bbox=[-72.5,40.5,-72,41]
    intersects=dict(type="Point", coordinates=[lon, lat]),
)
items = search.get_all_items()
items[0].assets  # looks correct, L1C products

stack = stackstac.stack(items, epsg=4326)

lowcloud = stack[stack["eo:cloud_cover"] < 20]
rgb = lowcloud.sel(band=["red", "green", "blue"])
monthly = rgb.resample(time="MS").median("time", keep_attrs=True)

x_utm, y_utm = pyproj.Proj(monthly.crs)(lon, lat)
buffer = 2000  # meters

aoi = monthly.loc[..., y_utm+buffer:y_utm-buffer, x_utm-buffer:x_utm+buffer]   # 2.65 GB Xarray, looks correct

aws_session = AWSSession(requester_pays=True)

with rio.Env(aws_session):
    with dask.diagnostics.ProgressBar():
        data = aoi.compute()

CPLE_OpenFailedError RuntimeError: Error opening 's3://sentinel-s2-l1c/tiles/13/S/DV/2023/7/20/0/B02.jp2': RasterioIOError("'/vsis3/sentinel-s2-l1c/tiles/13/S/DV/2023/7/20/0/B02.jp2' does not exist in the file system, and is not recognized as a supported dataset name.")

J6767 commented 1 year ago

Incidentally, the Boto version does find and download the products, but I want to subset them before downloading, using Stackstac.

import boto3

s3_client = boto3.Session().client('s3')
response = s3_client.get_object(Bucket='sentinel-s2-l1c',
                                Key='tiles/7/W/FR/2018/3/31/0/B01.jp2', 
                                RequestPayer='requester')
response_content = response['Body'].read()

with open('./B01.jp2', 'wb') as file:
    file.write(response_content)
J6767 commented 1 year ago

The error originated with my installation, works like a charm now. Thanks for your hard work on Stackstac!