#!/usr/bin/env python
# coding: utf-8

# ## Accessing MODIS snow data with the Planetary Computer STAC API
# 
# The planetary computer hosts two snow-related MODIS 6.1 products:
# 
# - Snow cover daily (10A1)
# - Snow cover 8-day (10A2)
# 
# For more information about the products themselves, check out the Data Pages at the [bottom of this document](#data-pages).

# ### Environment setup
# 
# This notebook works with or without an API key, but you will be given more permissive access to the data with an API key.
# The Planetary Computer Hub is pre-configured to use your API key.

# In[1]:


import odc.stac
import planetary_computer
import pystac_client
import rich.table


# ### Data access
# 
# The datasets hosted by the Planetary Computer are available from [Azure Blob Storage](https://docs.microsoft.com/en-us/azure/storage/blobs/). We'll use [pystac-client](https://pystac-client.readthedocs.io/) to search the Planetary Computer's [STAC API](https://planetarycomputer.microsoft.com/api/stac/v1/docs) for the subset of the data that we care about, and then we'll load the data directly from Azure Blob Storage. We'll specify a `modifier` so that we can access the data stored in the Planetary Computer's private Blob Storage Containers. See [Reading from the STAC API](https://planetarycomputer.microsoft.com/docs/quickstarts/reading-stac/) and [Using tokens for data access](https://planetarycomputer.microsoft.com/docs/concepts/sas/) for more.

# In[2]:


catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1",
    modifier=planetary_computer.sign_inplace,
)


# ### Query for available data
# 
# MODIS is a global dataset with a variety of products available within each larger category (vegetation, snow, fire, temperature, and reflectance).
# The [MODIS group](https://planetarycomputer.microsoft.com/dataset/group/modis) contains a complete listing of available collections.
# Each collection's id is in the format `modis-{product}-061`, where `product` is the MODIS product id.
# The `-061` suffix indicates that all of the MODIS collections are part of the [MODIS 6.1 update](https://atmosphere-imager.gsfc.nasa.gov/documentation/collection-61).
# 
# We'll look at the snow cover around the Mount of the Holy Cross in Colorado and how it progresses throughout the winter, using the daily snow product (10A1).

# In[3]:


longitude = -106.481687
latitude = 39.466829
mount_of_the_holy_cross = [longitude, latitude]
geometry = {
    "type": "Point",
    "coordinates": mount_of_the_holy_cross,
}
datetimes = [
    "2020-11",
    "2020-12",
    "2021-01",
    "2021-02",
    "2021-03",
    "2021-04",
    "2021-05",
    "2021-06",
]
items = dict()

for datetime in datetimes:
    print(f"Fetching {datetime}")
    search = catalog.search(
        collections=["modis-10A1-061"],
        intersects=geometry,
        datetime=datetime,
    )
    items[datetime] = search.get_all_items()[0]

print(items)


# ### Available assets
# 
# Each item has several available assets, including the original HDF file and a Cloud-optimized GeoTIFF of each subdataset.

# In[4]:


t = rich.table.Table("Key", "Title")
for key, asset in items["2020-11"].assets.items():
    t.add_row(key, asset.title)
t


# ### Loading the snow cover data
# 
# For this example, we'll visualize and compare the snow cover month to month.
# Let's grab each snow cover COG and load them into an xarray using [odc-stac](https://github.com/opendatacube/odc-stac).
# We'll also apply the scaling as defined by the `raster:bands` extension.
# The MODIS coordinate reference system is a [sinusoidal grid](https://modis-land.gsfc.nasa.gov/MODLAND_grid.html), which means that views in a naïve XY raster look skewed.
# For visualization purposes, we reproject to a [spherical Mercator projection](https://wiki.openstreetmap.org/wiki/EPSG:3857) for intuitive, north-up visualization.
# 
# The NDSI Snow Cover values are defined as:
# 
# ```
# 0–100: NDSI snow cover
# 200: missing data
# 201: no decision
# 211: night
# 237: inland water
# 239: ocean
# 250: cloud
# 254: detector saturated
# 255: fill
# ```
# 
# We want to mask out all numbers greater than 100.

# In[5]:


bbox = [longitude - 0.5, latitude - 0.5, longitude + 0.5, latitude + 0.5]
data = odc.stac.load(
    items.values(),
    crs="EPSG:3857",
    bbox=bbox,
    bands="NDSI_Snow_Cover",
    resolution=500,
)
data = data.where(data <= 100, drop=True)
data


# ### Displaying the data
# 
# Let's display the snow cover for each month.

# In[6]:


g = data["NDSI_Snow_Cover"].plot.imshow(
    col="time", vmin=0, vmax=100, col_wrap=4, size=4
)
months = data["NDSI_Snow_Cover"].time.to_pandas().dt.strftime("%B %Y")

for ax, month in zip(g.axes.flat, months):
    ax.set_title(month)


# You'll notice there's a lot of missing values due to masking, probably due to clouds.
# More sophisticated analysis would use the `NDSI_Snow_Cover_Basic_QA` and other bands to merge information from multiple scenes, or pick the best scenes for analysis.
# 
# ### Data pages
# 
# These pages include links to the user guides:
# 
# - MOD10A1: https://nsidc.org/data/MOD10A1/versions/61
# - MOD10A2: https://nsidc.org/data/MOD10A2/versions/61