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

# ## Accessing Landsat Collection 2 Level-1 and Level-2 data with the Planetary Computer STAC API
# 
# The Planetary Computer's [Landsat](https://landsat.gsfc.nasa.gov/) dataset represents a global archive of [Level-1](https://www.usgs.gov/landsat-missions/landsat-collection-2-level-1-data) and [Level-2](https://www.usgs.gov/landsat-missions/landsat-collection-2-level-2-science-products) data from from [Landsat Collection 2](https://www.usgs.gov/core-science-systems/nli/landsat/landsat-collection-2). The dataset is grouped into two STAC Collections:
# - `landsat-c2-l2` for Level-2 data
# - `landsat-c2-l1` for Level-1 data
# 
# This notebook demonstrates the use of the Planetary Computer STAC API to query for Landsat data. We will start with an example using Level-2 data and follow with a similar example using Level-1 data.
# 
# ### 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](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 `pc.settings.set_subscription_key(<YOUR API Key>)`.

# In[1]:


import pystac_client
import planetary_computer
import odc.stac
import matplotlib.pyplot as plt

from pystac.extensions.eo import EOExtension as eo


# ### 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,
)


# ### Choose an area and time of interest
# 
# This area is in Redmond, WA, USA, near Microsoft's main campus. We'll search all of 2021, limiting the results to those less than 10% cloudy.

# In[3]:


bbox_of_interest = [-122.2751, 47.5469, -121.9613, 47.7458]
time_of_interest = "2021-01-01/2021-12-31"


# In[4]:


search = catalog.search(
    collections=["landsat-c2-l2"],
    bbox=bbox_of_interest,
    datetime=time_of_interest,
    query={"eo:cloud_cover": {"lt": 10}},
)

items = search.item_collection()
print(f"Returned {len(items)} Items")


# Let's find the least cloudy of the bunch.

# In[5]:


selected_item = min(items, key=lambda item: eo.ext(item).cloud_cover)

print(
    f"Choosing {selected_item.id} from {selected_item.datetime.date()}"
    + f" with {selected_item.properties['eo:cloud_cover']}% cloud cover"
)


# ### Available assets
# 
# In addition to numerous metadata assets, each Electro-Optical (EO) band is a separate asset.

# In[6]:


max_key_length = len(max(selected_item.assets, key=len))
for key, asset in selected_item.assets.items():
    print(f"{key.rjust(max_key_length)}: {asset.title}")


# ### Render a natural color image of the AOI
# 
# We'll start by loading the red, green, and blue bands for our area of interest into an xarray dataset using [`odc-stac`](https://pypi.org/project/odc-stac/). We will also load the nir08 band for use in computing an NDVI value in a later example. Note that we pass the `sign` function from the [planetary-computer](https://github.com/microsoft/planetary-computer-sdk-for-python) package so that `odc-stac` can supply the required [Shared Access Signature](https://docs.microsoft.com/en-us/azure/storage/common/storage-sas-overview) tokens that are necessary to download the asset data.

# In[7]:


bands_of_interest = ["nir08", "red", "green", "blue", "qa_pixel", "lwir11"]
data = odc.stac.stac_load(
    [selected_item], bands=bands_of_interest, bbox=bbox_of_interest
).isel(time=0)
data


# Now we'll convert our xarray Dataset to a DataArray and plot the RGB image. We set `robust=True` in `imshow` to avoid manual computation of the color limits (vmin and vmax) that is necessary for data not scaled to between 0-1 while also eliminating extreme values that can cause a washed out image.

# In[8]:


fig, ax = plt.subplots(figsize=(10, 10))

data[["red", "green", "blue"]].to_array().plot.imshow(robust=True, ax=ax)
ax.set_title("Natural Color, Redmond, WA");


# ### Render an NDVI image of the AOI
# 
# Landsat has several bands, and with them we can go beyond rendering natural color imagery; for example, the following code computes a [Normalized Difference Vegetation Index (NDVI)](https://en.wikipedia.org/wiki/Normalized_difference_vegetation_index) using the near-infrared and red bands. Note that we convert the red and near infrared bands to a data type that can contain negative values; if this is not done, negative NDVI values will be incorrectly stored.

# In[9]:


red = data["red"].astype("float")
nir = data["nir08"].astype("float")
ndvi = (nir - red) / (nir + red)

fig, ax = plt.subplots(figsize=(14, 10))
ndvi.plot.imshow(ax=ax, cmap="viridis")
ax.set_title("NDVI, Redmond, WA");


# ### Selecting specific platforms
# 
# Landsat Collection 2 Level-2 includes assets from several different Platforms.

# In[10]:


catalog.get_collection("landsat-c2-l2").summaries.to_dict()["platform"]


# You might want to limit your search to a specific platform or platforms, to avoid the [Landsat 7 Scan Line Corrector failure](https://www.usgs.gov/landsat-missions/landsat-7), for example. Use the `"platform": {"in": ... }` query to select specific platforms.

# In[11]:


search = catalog.search(
    collections=["landsat-c2-l2"],
    bbox=bbox_of_interest,
    datetime=time_of_interest,
    query={
        "eo:cloud_cover": {"lt": 10},
        "platform": {"in": ["landsat-8", "landsat-9"]},
    },
)
items = search.get_all_items()


# ### Rescaling Temperature Data
# 
# Landsat Collection 2 Level 2 includes measures of [Surface Temperature](https://www.usgs.gov/landsat-missions/landsat-surface-temperature). We'll look at the surface temperature band `TIRS_B10` under the key `lwir11`. The raw values are rescaled, so you should scale and offset the data before interpreting it. Use the metadata in the asset's `raster_bands` to find the scale and offset values:

# In[12]:


band_info = selected_item.assets["lwir11"].extra_fields["raster:bands"][0]
band_info


# To go from the raw values to Kelvin, we apply those values:

# In[13]:


temperature = data["lwir11"].astype(float)
temperature *= band_info["scale"]
temperature += band_info["offset"]
temperature[:5, :5]


# To convert from Kelvin to degrees, subtract 273.15.

# In[14]:


celsius = temperature - 273.15
celsius.plot(cmap="magma", size=10);


# ### Landsat Collection 2 Level-1 Data
# 
# Thus far we have worked with Landsat Collection 2 Level-2 data, which is processed to a consistent set of surface reflectance and surface temperature [science products](https://www.usgs.gov/landsat-missions/landsat-science-products). 
# 
# The Planetary Computer also hosts Collection 2 Level-1 data (top of atmosphere values) acquired by the [Multispectral Scanner System](https://landsat.gsfc.nasa.gov/multispectral-scanner-system/) (MSS) onboard Landsat 1 through 5. These data do not include a blue band, so a natural color image is not possible. We will plot a color infrared image from the nir08, red, and green bands instead.
# 
# As before, we use pystac-client to search over the `landsat-c2-l1` collection. We'll use the same area of interest and return only those with less than 10% cloud cover.

# In[15]:


search = catalog.search(
    collections=["landsat-c2-l1"],
    bbox=bbox_of_interest,
    query={"eo:cloud_cover": {"lt": 10}},
)

items = search.item_collection()
print(f"Returned {len(items)} Items")


# Choose the least cloudy Item.

# In[16]:


selected_item = min(items, key=lambda item: eo.ext(item).cloud_cover)

print(
    f"Choosing {selected_item.id} from {selected_item.datetime.date()}"
    + f" with {selected_item.properties['eo:cloud_cover']}% cloud cover"
)


# ### Available assets
# 
# MSS data has four EO bands assets and a number of metadata files and bands.

# In[17]:


max_key_length = len(max(selected_item.assets, key=len))
for key, asset in selected_item.assets.items():
    print(f"{key.rjust(max_key_length)}: {asset.title}")


# We'll use `odc-stac` to load only the EO bands and area we are interested in.

# In[18]:


bands_of_interest = ["nir08", "red", "green"]
data = odc.stac.stac_load(
    [selected_item], bands=bands_of_interest, bbox=bbox_of_interest
).isel(time=0)

cir = data.to_array()

fig, ax = plt.subplots(figsize=(10, 10))
cir.plot.imshow(robust=True, ax=ax)
ax.set_title("Color Infrared, Redmond, WA");