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

# ## Using ALOS PALSAR and Forest/Non-Forest Annual Mosaics with the Planetary Computer STAC API

# ALOS PALSAR Annual Mosaic consists of radar backscatter in the HH and HV with additional bands for ancillary information about incidince angle, masking, and date of acquisition (per pixel).
# 
# For this example we'll find a 1x1 degree data tile with a diversity of values, and plot the HH, HV and matching Forest Classification.

# In[1]:


import numpy as np
import planetary_computer
import stackstac
import pystac_client
import warnings

warnings.filterwarnings("ignore", "The argument 'infer_datetime_format'", UserWarning)


# ### Polarization Returns
# 
# We'll use the Planetary Computer's STAC API to find some scenes matching our area of interest. See [reading STAC](https://planetarycomputer.microsoft.com/docs/quickstarts/reading-stac/) for more.

# In[2]:


bbox = [9.4, 0, 9.5, 1]
client = pystac_client.Client.open(
    "https://planetarycomputer.microsoft.com/api/stac/v1/",
    modifier=planetary_computer.sign_inplace,
)
search = client.search(
    collections=["alos-palsar-mosaic"],
    bbox=bbox,
)
items = search.item_collection()
print(f"Returned {len(items)} items")


# Each of these items has a handful of assets.

# In[3]:


item = items[0]
list(item.assets)


# We'll select the first item and load it into xarray using stackstac. Note that the underlying data have a `uint16` dtype, so we specify that when creating the `DataArray`.

# In[4]:


data = stackstac.stack(
    items[0],
    assets=["HH", "HV"],  # Both Polarizations
    dtype="uint16",
    fill_value=0,
).squeeze()
data


# Most of the values are near zero, but there's a long tail.

# In[5]:


data.sel(band="HH").plot.hist(bins=100, size=5);


# In[6]:


img = data.sel(band="HV").plot(vmin=0, vmax=14000, size=8)
img.axes.set_axis_off();


# ### Forest Classification
# 
# ALOS Forest/Non-Forest Classification is derived from the ALOS PALSAR Annual Mosaic, and classifies the pixels to detect forest cover. We'll search for Forest / Non-Forest items covering the same area.

# In[7]:


import shapely.geometry

s = shapely.geometry.shape(item.geometry)
intersects = shapely.geometry.mapping(s.centroid)

search = client.search(collections=["alos-fnf-mosaic"], intersects=intersects)
fnf_items = search.get_all_items()
print(f"Returned {len(fnf_items)} items")


# The primary asset is under the key `C`.

# In[8]:


fnf_item = fnf_items[0]

fnf_tile = stackstac.stack(
    fnf_item, assets=["C"], dtype="uint8", fill_value=0
).squeeze()
fnf_tile


# We can get the class name from the `classification:classes` field on the Collection's `item_assets`.

# In[9]:


collection = client.get_collection("alos-fnf-mosaic")

classes = collection.extra_fields["item_assets"]["C"]["classification:classes"]
labels = [x["description"] for x in classes]


# In[10]:


import matplotlib.colors

cmap = matplotlib.colors.ListedColormap(
    [
        (0, 0, 0, 0),  # nodata
        (0, 0.7, 0, 1),  # forest, > 90% canopy
        (0.51, 0.94, 0.38, 1),  # forest 10-90% canopy
        (1, 1, 0.6, 1),  # non-forest
        (0.0, 0.0, 1, 1),  # water
    ]
)

# Plot the classification map
p = fnf_tile.plot(cmap=cmap, vmin=0, vmax=4, size=8)
ticks = np.arange(0.4, 4, 0.8)
p.colorbar.set_ticks(ticks, labels=labels)
p.colorbar.set_label("Class")