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

# ## Accessing Sentinel-3 OLCI data with the Planetary Computer STAC API
# 
# The [Sentinel 3 OLCI instrument](https://sentinel.esa.int/web/sentinel/technical-guides/sentinel-3-olci) provides radiance measurements of the Earth's surface in the visible and near infra-red spectral domain.
# These measurements are processed into two Level 2 products, one for the ocean and one for the land.
# Each product has its own STAC collection in the Planetary Computer:
# 
# - Land: `sentinel-3-olci-lfr-l2-netcdf`
# - Ocean: `sentinel-3-olci-wfr-l2-netcdf`
# 
# This notebook demonstrates accessing and visualizing data from both collections.
# 
# ### Data Access
# 
# This notebook works with or without an API key, but you will be given more permissive access to the data with an API key. If you are using the [Planetary Computer Hub](https://planetarycomputer.microsoft.com/compute) to run this notebook, then your API key is automatically set to the environment variable `PC_SDK_SUBSCRIPTION_KEY` for you when your server is started. Otherwise, you can view your keys by signing in to the [developer portal](https://planetarycomputer.developer.azure-api.net/). The API key may be manually set via the environment variable `PC_SDK_SUBSCRIPTION_KEY` or the following code:
# 
# ```python
# import planetary_computer
# planetary_computer.settings.set_subscription_key(<YOUR API Key>)
# ```
# 
# The datasets hosted by the Planetary Computer are available in [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.
# 
# 

# ### Water products
# 
# The collection description tells us more about the water products.

# In[1]:


import planetary_computer
import pystac_client
from IPython.display import display, Markdown

catalog = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1",
    modifier=planetary_computer.sign_inplace,
)
collection = catalog.get_collection("sentinel-3-olci-wfr-l2-netcdf")
display(Markdown(f"### {collection.id}\n\n{collection.description}"))


# ### Define the area of interest and search the water collection
# 
# We'll search for items over the coordinates `[-2.79, 44.14]`.

# In[2]:


import xarray as xr
import fsspec

search = catalog.search(
    collections=["sentinel-3-olci-wfr-l2-netcdf"],
    intersects={"type": "Point", "coordinates": [-2.79, 44.14]},
)
item = next(search.items())


# ### Available Assets and Metadata
# 
# Each item includes a handful of assets linking to NetCDF files with the data or additional metadata files.

# In[3]:


import rich.table

t = rich.table.Table("Key", "Value")
for key, asset in item.assets.items():
    t.add_row(key, asset.description)

t


# ### Reading data
# 
# We can use xarray to read each NetCDF file directly from Blob Storage.

# In[4]:


keys = [
    "iwv",
    "par",
    "trsp",
    "w-aer",
    "chl-nn",
    "iop-nn",
    "tsm-nn",
    "chl-oc4me",
    "oa01-reflectance",
]
datasets = [xr.open_dataset(fsspec.open(item.assets[k].href).open()) for k in keys]

ds = xr.combine_by_coords(datasets, join="exact", combine_attrs="drop_conflicts")
ds


# We'll plot the integrated water vapor variable.

# In[5]:


ds.IWV.coarsen({"rows": 10, "columns": 10}, boundary="trim").mean().plot();


# ### Geolocating the data
# 
# The geospatial information in this dataset is distributed as a separate NetCDF file, cataloged under the `geoCoordinates` key. That contains a dataset with `latitude` and `longitude` arrays, each of which is the same shape as the data variables and gives the latitude and longitude for each pixel in data variable.
# 
# We'll reshape the data to a (long-form) DataFrame with a single row for each pixel. We'll then make a scatter plot, using the longitude and latitude as the `x` and `y` coordintes.

# In[6]:


import pandas as pd
import datashader
import colorcet

geo = xr.open_dataset(fsspec.open(item.assets["geo-coordinates"].href).open()).load()
a01 = ds.Oa01_reflectance.load()


df = pd.DataFrame(
    {
        "longitude": geo.longitude.data.ravel(),
        "latitude": geo.latitude.data.ravel(),
        "value": a01.data.ravel(),
    }
)


# To avoid overplotting the data, we'll use [datashader](https://datashader.org/).

# In[7]:


cvs = datashader.Canvas(plot_width=600, plot_height=600)
agg = cvs.points(df, "latitude", "longitude", agg=datashader.reductions.mean("value"))
img = datashader.tf.shade(agg, cmap=colorcet.CET_CBD2)
img


# ### Land data
# 
# Let's do the same process, but for the land product.

# In[8]:


collection = catalog.get_collection("sentinel-3-olci-lfr-l2-netcdf")
display(Markdown(f"### {collection.id}\n\n{collection.description}"))


# In[9]:


search = catalog.search(
    collections=["sentinel-3-olci-lfr-l2-netcdf"],
    intersects={"type": "Point", "coordinates": [-2.79, 44.14]},
)
item = next(search.items())
t = rich.table.Table("Key", "Value")
for key, asset in item.assets.items():
    t.add_row(key, asset.description)

t


# ### Plot FAPAR
# 
# The green instantaneous Fraction of Absorbed Photosynthetically Active Radiation (FAPAR) product uses a [vegetation index algorithm](https://sentinel.esa.int/web/sentinel/technical-guides/sentinel-3-olci/level-2/olci-global-vegetation-index) to provide an estimate of the FARAR in the plant canopy.

# In[10]:


dataset = xr.open_dataset(fsspec.open(item.assets["gifapar"].href).open())
dataset


# In[11]:


geo = xr.open_dataset(fsspec.open(item.assets["geo-coordinates"].href).open()).load()
fapar = dataset.GIFAPAR.load()

df = pd.DataFrame(
    {
        "longitude": geo.longitude.data.ravel(),
        "latitude": geo.latitude.data.ravel(),
        "value": fapar.data.ravel(),
    }
)
cvs = datashader.Canvas(plot_width=600, plot_height=600)
agg = cvs.points(df, "latitude", "longitude", agg=datashader.reductions.mean("value"))
img = datashader.tf.shade(agg, cmap=colorcet.CET_CBD2)
img