GET NDVI from Element84 using pystac

ndvi.py

from pystac_client import Client
from odc.stac import load
import odc.geo
# Load the STAC catalog

client = Client.open("https://earth-search.aws.element84.com/v1") 

# load the STAC collection
collection = client.get_collection("sentinel-2-c1-l2a")

# Geometry for the area of interest
geometry = {
        "coordinates": [
          [
            [
              73.79383362280944,
              19.94813852484772
            ],
            [
              73.79383362280944,
              19.94328417001347
            ],
            [
              73.79808489158432,
              19.94328417001347
            ],
            [
              73.79808489158432,
              19.94813852484772
            ],
            [
              73.79383362280944,
              19.94813852484772
            ]
          ]
        ],
        "type": "Polygon"
      }

# Add more filtering parameters
filtering_params = {
    "eo:cloud_cover": {"lt": 10},  # Less than 20
}
# Load the STAC catalog for a particular date 
date_YYMMDD = "2024-01-01/2024-01-15"
# run pystac client search to see available dataset
search = client.search(
    collections=[collection], intersects=geometry, datetime=date_YYMMDD, query=filtering_params
)
# spit out data as GeoJSON dictionary
print(search.item_collection_as_dict())

for item in search.items():
    print(item.id)
    
data = load(search.items() ,geopolygon=geometry,groupby="solar_day", chunks={})
# Applying scale and offset

def getstuff(band,prop):
  return search.item_collection_as_dict()['features'][0]['assets'][band]['raster:bands'][0][prop]
data["red"] = data["red"] * getstuff('red', 'scale') + getstuff('red', 'offset')
data["nir"] = data["nir"] * getstuff('nir', 'scale') + getstuff('nir', 'offset')

ndvi = (data.nir - data.red) / (data.nir + data.red)
ndvi_1 = ndvi.isel(time=1)
odc.geo.xr.write_cog(ndvi_1,fname='ndvi_1_scaled.tiff',  overwrite=True)

req.txt

pystac-client
odc-stac
rioxarray