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

# ## Accessing Sentinel-1 GRD data with the Planetary Computer STAC API
# 
# [Sentinel-1 Ground Range Detected (GRD)](https://planetarycomputer.microsoft.com/dataset/sentinel-1-grd) Level-1 products produced by the European Space Agency.
# 
# ### Environment setup
# 
# Running this notebook requires an API key.
# 
# * The [Planetary Computer Hub](https://planetarycomputer.microsoft.com/compute) is pre-configured to use your API key.
# * To use your API key locally, set the environment variable `PC_SDK_SUBSCRIPTION_KEY` or use `planetary_computer.settings.set_subscription_key(<YOUR API Key>)`
# 
# See [when an account is needed](https://planetarycomputer.microsoft.com/docs/concepts/sas/#when-an-account-is-needed) for more, and [request an account](http://planetarycomputer.microsoft.com/account/request) if needed.

# In[1]:


import ipyleaflet
import matplotlib.pyplot as plt
import numpy as np
import pystac
import pystac_client
import planetary_computer
import requests
import rich.table
import xarray as xr
import odc.stac
import hvplot.xarray

from IPython.display import Image
from pystac_client.stac_api_io import StacApiIO
from pystac_client.client import Client


# Because NSF uses a self-signed certificate, we **can't** use the standard method for authenticating with the Planetary Computer:
# ```
# catalog = pystac_client.Client.open(
#     "https://planetarycomputer.microsoft.com/api/stac/v1",
#     modifier=planetary_computer.sign_inplace,
# )
# ```
# instead, we need to point to our custom certificate thusly:

# In[2]:


stac_api_io = StacApiIO()
stac_api_io.session.verify = "/home/shared/PA-RootCA-Cert-2023-Pub.pem"

catalog = Client.open("https://planetarycomputer.microsoft.com/api/stac/v1", stac_io=stac_api_io)


# ### 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.

# ### Choose an area and time of interest
# 
# We'll search for assets acquired in a region of the Arctic in July 2023. 

# In[3]:


collection = 'sentinel-1-grd'


# In[4]:


#bbox = [-80.11, 8.71, -79.24, 9.38]
bbox = [1.5, 78.5, 3.5, 79.1]
search = catalog.search(
    collections=[collection], bbox=bbox, datetime="2023-07-03/2023-07-11"
)
items = search.item_collection()
print(f"Found {len(items)} items")
item = items[0]


# In[5]:


for item in items:
    print(item)


# ### Inspect the STAC metadata
# 
# The STAC metadata includes many useful pieces of metadata, including metadata from the [SAR](https://github.com/stac-extensions/sar) and [Satellite](https://github.com/stac-extensions/sat) extensions.

# In[6]:


table = rich.table.Table("key", "value")
for k, v in sorted(item.properties.items()):
    table.add_row(k, str(v))

table


# The data assets on every Sentinel-1 RTC item will be some combination of `hh`, `hv`, `vh`, and `vv`. These represent the terrain-corrected gamma nought values of a signal transmitted in one polarization ("h" or "v") and received in another ("h" or "v"). The `sar:polarizations` field indicates which assets are available.

# In[7]:


item.properties["sar:polarizations"]


# ### Visualize the assets
# 
# Next, we'll load the `vv` data into [xarray](https://xarray.pydata.org/) and plot the results. We'll use [Dask](http://dask.org/) to load the data in parallel. We're working with a small amount of data so we'll use a single machine. For larger datasets, see [Scaling with Dask](https://planetarycomputer.microsoft.com/docs/quickstarts/scale-with-dask/).

# In[8]:


from distributed import Client

client = Client(processes=False)
print(client.dashboard_link)


# #### Load data with `odc.stac`

# In[9]:


sas_token = stac_api_io.session.get(f"https://planetarycomputer.microsoft.com/api/sas/v1/token/{collection}").json()["token"


# We need to append the `sas_token` to each of the `href` URLS so that stackstack or odc.stac will work:

# In[10]:


for i in range(len(items)):
    for asset in items[i].assets.keys():
        items[i].assets[asset].href = f"{items[i].assets[asset].href}?{sas_token}"


# In[12]:


bands_of_interest = ['hh', 'hv']

da = odc.stac.stac_load(
    items,
    bands=bands_of_interest,
    bbox=bbox,
    crs="EPSG:3857",
    resolution=100,
    chunks={},  # <-- use Dask
).to_array(dim='band')
da


# In[13]:


da.sel(band='hh').hvplot.image(x='x', y='y', crs=3857, rasterize=True, tiles='OSM')


# ### Can't get stackstac to work

# In[14]:


import stackstac

da = stackstac.stack(items, bounds_latlon=bbox, epsg=32630, resolution=100)
da


# In[15]:


da.sel(band="hh")[0]


# We'll select the `vv` band for the first timestep found by our search.

# In[16]:


hh = da.sel(band="hh")[0].compute()


# The distribution of the raw values is quite skewed:

# In[ ]:


hh.plot.hist(bins=30);


# So the values are typically transformed before visualization:

# In[ ]:


fig, ax = plt.subplots(figsize=(6, 4))


def db_scale(x):
    return 10 * np.log10(x)


db_scale(vv).plot.hist(bins=50, ax=ax)
ax.set(title="Distribution of pixel values (dB scale)", xlabel="Pixel values");


# In[ ]:


img = (
    db_scale(vv)
    .coarsen(x=4, y=4, boundary="trim")
    .max()
    .plot.imshow(cmap="bone", size=12, aspect=1.05, add_colorbar=False)
)
img.axes.set_axis_off();


# ### The effect of terrain correction
# 
# In this section, we compare Sentinel-1 GRD to Sentinel-1 RTC to see the effect of terrain correction.
# 
# Every Sentinel-1-RTC item is derived from a [Sentinel-1-GRD](https://planetarycomputer.microsoft.com/dataset/sentinel-1-grd) item. You can follow the `derived_from` link to get back to the original GRD item.

# In[ ]:


rtc_item = catalog.get_collection("sentinel-1-rtc").get_item(
    "S1A_IW_GRDH_1SDV_20220518T054334_20220518T054359_043261_052A9D_rtc"
)
grd_item = pystac.read_file(rtc_item.get_single_link("derived_from").target)


# Next, we'll use the `tilejson` asset, which uses the Planetary Computer's [Data API](https://planetarycomputer.microsoft.com/api/data/v1/) to serve xyz tiles for a STAC item.

# In[ ]:


grd_tiles = requests.get(grd_item.assets["tilejson"].href).json()["tiles"][0]
rtc_tiles = requests.get(rtc_item.assets["tilejson"].href).json()["tiles"][0]


# With these URLs, we can build an interactive map using [ipyleaflet](https://ipyleaflet.readthedocs.io/en/latest/index.html). Adjust the slider to visualize either GRD (to the left) or RTC (to the right).

# In[ ]:


center = [47.05, 7.10]
m = ipyleaflet.Map(
    center=center,
    zoom=14,
    controls=[ipyleaflet.FullScreenControl()],
)
grd_layer = ipyleaflet.TileLayer(url=grd_tiles)
rtc_layer = ipyleaflet.TileLayer(url=rtc_tiles)

control = ipyleaflet.SplitMapControl(left_layer=grd_layer, right_layer=rtc_layer)
m.add_control(control)
m.scroll_wheel_zoom = True
m


# Notice that points seem to "jump" between the GRD and RTC. The RTC values are corrected to align with where they're actually at on the Earth.
# 
# For more background on terrain correction, and for an introduction to the [sarsen](https://github.com/bopen/sarsen) package which enables customizable RTCs, see [Sentinel-1 Customizable Radiometric Terrain Correction](https://planetarycomputer.microsoft.com/docs/tutorials/customizable-rtc-sentinel1/).