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 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")

Returned 21 items

Each of these items has a handful of assets.

In [3]:

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

Out[3]:

['HH',
 'HV',
 'date',
 'mask',
 'linci',
 'metadata',
 'tilejson',
 'rendered_preview']

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

Out[4]:

<xarray.DataArray 'stackstac-1a4421f08bc5ca33d576ce09d60ce330' (band: 2,
                                                                y: 4501, x: 4500)>
dask.array<getitem, shape=(2, 4501, 4500), dtype=uint16, chunksize=(1, 1024, 1024), chunktype=numpy.ndarray>
Coordinates: (12/27)
    time                            datetime64[ns] 2021-01-01
    id                              <U21 'N02E009_21_F02DAR_MOS'
  * band                            (band) <U2 'HH' 'HV'
  * x                               (x) float64 9.0 9.0 9.0 ... 9.999 10.0 10.0
  * y                               (y) float64 2.0 2.0 2.0 ... 1.001 1.0 1.0
    proj:epsg                       int64 4326
    ...                              ...
    start_datetime                  <U20 '2021-01-01T00:00:00Z'
    platform                        <U6 'ALOS-2'
    gsd                             int64 25
    proj:bbox                       object {9.0, 10.0, 2.0000000000000004, 1....
    raster:bands                    object {'nodata': 1, 'data_type': 'uint16'}
    epsg                            int64 4326
Attributes:
    spec:        RasterSpec(epsg=4326, bounds=(9.0, 1.0, 10.0, 2.000222222222...
    crs:         epsg:4326
    transform:   | 0.00, 0.00, 9.00|\n| 0.00,-0.00, 2.00|\n| 0.00, 0.00, 1.00|
    resolution:  0.00022222222222222223

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")

Returned 6 items

/srv/conda/envs/notebook/lib/python3.10/site-packages/pystac_client/item_search.py:841: FutureWarning: get_all_items() is deprecated, use item_collection() instead.
  warnings.warn(

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

Out[8]:

<xarray.DataArray 'stackstac-3882cecac9b4ee67fd63e1c7d3c237f3' (y: 4500, x: 4500)>
dask.array<getitem, shape=(4500, 4500), dtype=uint8, chunksize=(1024, 1024), chunktype=numpy.ndarray>
Coordinates: (12/18)
    time            datetime64[ns] 2020-01-01
    id              <U14 'N02E009_20_FNF'
    band            <U1 'C'
  * x               (x) float64 9.0 9.0 9.0 9.001 ... 9.999 9.999 10.0 10.0
  * y               (y) float64 2.0 2.0 2.0 1.999 1.999 ... 1.001 1.001 1.0 1.0
    proj:shape      object {4500}
    ...              ...
    instruments     <U8 'PALSAR-2'
    title           <U1 'C'
    gsd             int64 25
    proj:bbox       object {9.0, 10.0, 2.0, 1.0}
    raster:bands    object {'nodata': 0, 'data_type': 'uint8'}
    epsg            int64 4326
Attributes:
    spec:        RasterSpec(epsg=4326, bounds=(9.0, 1.0, 10.0, 2.0), resoluti...
    crs:         epsg:4326
    transform:   | 0.00, 0.00, 9.00|\n| 0.00,-0.00, 2.00|\n| 0.00, 0.00, 1.00|
    resolution:  0.00022222222222222223

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")