import geopandas as gpd
import planetary_computer
import pystac_client
import odc.stac
import numpy as np

bbox = [105.5, 9.2, 106.3, 10.0]
time_range = "2023-01-01/2023-04-30"
catalog = pystac_client.Client.open("https://planetarycomputer.microsoft.com/api/stac/v1", modifier=planetary_computer.sign_inplace)
items = list(catalog.search(collections=["sentinel-2-l2a"], bbox=bbox, datetime=time_range, query={"eo:cloud_cover": {"lt": 30}}).items())
items = [planetary_computer.sign(item) for item in items]

x = 561609
y = 1024183
ds = odc.stac.load(items, bands=["B02"], x=(x-80, x+80), y=(y-80, y+80), crs="EPSG:32648", resolution=10, patch_url=planetary_computer.sign, fail_on_error=False).compute()
print("Original time size:", len(ds.time))
ds2 = ds.dropna(dim="time", how="all")
print("After dropna time size:", len(ds2.time))
print("B02 mean:", np.nanmean(ds2["B02"].values))
print("B02 non-nan count:", np.sum(~np.isnan(ds2["B02"].values)))