Mapping Cyanobacteria with PACE data¶

Install packages¶

Uncomment the following cell to install the HyperCoast package.

In [1]:

# %pip install -U "hypercoast[extra]"

# %pip install -U "hypercoast[extra]"

Import libraries¶

In [2]:

import earthaccess
import hypercoast

import earthaccess import hypercoast

Download PACE data¶

To download and access the PACE AOP data, you will need to create an Earthdata login. You can register for an account at urs.earthdata.nasa.gov. Once you have an account, run the following cell and enter your NASA Earthdata login credentials.

In [3]:

earthaccess.login(persist=True)

earthaccess.login(persist=True)

Out[3]:

<earthaccess.auth.Auth at 0x7f8b384f2c90>

Search for PACE AOP data:

In [4]:

results = hypercoast.search_pace(
    bounding_box=(-83, 25, -81, 28),
    temporal=("2024-07-30", "2024-08-15"),
    short_name="PACE_OCI_L2_AOP_NRT",
    count=1,
)

results = hypercoast.search_pace( bounding_box=(-83, 25, -81, 28), temporal=("2024-07-30", "2024-08-15"), short_name="PACE_OCI_L2_AOP_NRT", count=1, )

Download PACE AOP data:

In [5]:

hypercoast.download_pace(results[:1], out_dir="data")

hypercoast.download_pace(results[:1], out_dir="data")

Read PACE data¶

Read PACE AOP data as an xarray.Dataset:

In [6]:

filepath = "data/PACE_OCI.20240730T181157.L2.OC_AOP.V2_0.NRT.nc"
dataset = hypercoast.read_pace(filepath)
# dataset

filepath = "data/PACE_OCI.20240730T181157.L2.OC_AOP.V2_0.NRT.nc" dataset = hypercoast.read_pace(filepath) # dataset

Compute band ratios¶

In [7]:

da = dataset["Rrs"]
data = (
    (da.sel(wavelength=650) > da.sel(wavelength=620))
    & (da.sel(wavelength=701) > da.sel(wavelength=681))
    & (da.sel(wavelength=701) > da.sel(wavelength=450))
)
# data

da = dataset["Rrs"] data = ( (da.sel(wavelength=650) > da.sel(wavelength=620)) & (da.sel(wavelength=701) > da.sel(wavelength=681)) & (da.sel(wavelength=701) > da.sel(wavelength=450)) ) # data

The spectra of cyanobacteria bloom:¶

Visualize the selected region based on band ratios¶

In [8]:

import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature

# Create a plot
fig, ax = plt.subplots(figsize=(12, 6), subplot_kw={"projection": ccrs.PlateCarree()})

# ax.set_extent([-93, -87, 28, 32], crs=ccrs.PlateCarree())

# Plot the data
data.plot(
    ax=ax,
    transform=ccrs.PlateCarree(),
    cmap="coolwarm",
    cbar_kwargs={"label": "Cyano"},
)

# Add coastlines
ax.coastlines()

# Add state boundaries
states_provinces = cfeature.NaturalEarthFeature(
    category="cultural",
    name="admin_1_states_provinces_lines",
    scale="50m",
    facecolor="none",
)

ax.add_feature(states_provinces, edgecolor="gray")

# Optionally, add gridlines, labels, etc.
ax.gridlines(draw_labels=True)

plt.show()

import xarray as xr import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature # Create a plot fig, ax = plt.subplots(figsize=(12, 6), subplot_kw={"projection": ccrs.PlateCarree()}) # ax.set_extent([-93, -87, 28, 32], crs=ccrs.PlateCarree()) # Plot the data data.plot( ax=ax, transform=ccrs.PlateCarree(), cmap="coolwarm", cbar_kwargs={"label": "Cyano"}, ) # Add coastlines ax.coastlines() # Add state boundaries states_provinces = cfeature.NaturalEarthFeature( category="cultural", name="admin_1_states_provinces_lines", scale="50m", facecolor="none", ) ax.add_feature(states_provinces, edgecolor="gray") # Optionally, add gridlines, labels, etc. ax.gridlines(draw_labels=True) plt.show()

Cyanobacteria and Spectral Angle Mapper¶

Spectral Angle Mapper: Spectral similarity Input: library of Cyanobacteria bloom Rrs spectra with Chla at different levels

Spectral Mixture Analysis: unmix different cyanobacteria species based on spectral difference.

K-means applied to the whole image¶

In [9]:

import numpy as np
from sklearn.cluster import KMeans

import numpy as np from sklearn.cluster import KMeans

In [10]:

# Get the shape of the DataArray
print("Shape of da:", da.shape)

# Get the shape of the DataArray print("Shape of da:", da.shape)

Shape of da: (1710, 1272, 184)

In [11]:

# Get the dimension names of the DataArray
print("Dimension names:", da.dims)

# Get the dimension names of the DataArray print("Dimension names:", da.dims)

Dimension names: ('latitude', 'longitude', 'wavelength')

In [12]:

# Get the size of each dimension
print("Size of each dimension:", da.sizes)

# Get the size of each dimension print("Size of each dimension:", da.sizes)

Size of each dimension: Frozen({'latitude': 1710, 'longitude': 1272, 'wavelength': 184})

In [13]:

reshaped_data = da.values.reshape(-1, da.shape[-1])

reshaped_data = da.values.reshape(-1, da.shape[-1])

In [14]:

reshaped_data_no_nan = reshaped_data[~np.isnan(reshaped_data).any(axis=1)]

reshaped_data_no_nan = reshaped_data[~np.isnan(reshaped_data).any(axis=1)]

In [15]:

# Apply K-means clustering to classify into 5-6 water types.
n_clusters = 6
kmeans = KMeans(n_clusters=n_clusters, random_state=0)
kmeans.fit(reshaped_data_no_nan)

# Apply K-means clustering to classify into 5-6 water types. n_clusters = 6 kmeans = KMeans(n_clusters=n_clusters, random_state=0) kmeans.fit(reshaped_data_no_nan)

Out[15]:

KMeans(n_clusters=6, random_state=0)

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

In [16]:

# Initialize an array for cluster labels with NaN
labels = np.full(reshaped_data.shape[0], np.nan)

# Initialize an array for cluster labels with NaN labels = np.full(reshaped_data.shape[0], np.nan)

In [17]:

# Assign the computed cluster labels to the non-NaN positions
labels[~np.isnan(reshaped_data).any(axis=1)] = kmeans.labels_

# Assign the computed cluster labels to the non-NaN positions labels[~np.isnan(reshaped_data).any(axis=1)] = kmeans.labels_

In [18]:

# Reshape the labels back to the original spatial dimensions
cluster_labels = labels.reshape(da.shape[:-1])

# Reshape the labels back to the original spatial dimensions cluster_labels = labels.reshape(da.shape[:-1])

In [19]:

import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import cartopy.crs as ccrs
import cartopy.feature as cfeature

import matplotlib.pyplot as plt import matplotlib.colors as mcolors import cartopy.crs as ccrs import cartopy.feature as cfeature

In [20]:

# Assume 'cluster_labels' contains the K-means classification results reshaped to the original spatial dimensions
# Also assume that 'da' has the original latitude and longitude data
latitudes = da.coords["latitude"].values
longitudes = da.coords["longitude"].values

# Create a custom discrete color map for K-means clusters
cmap = mcolors.ListedColormap(
    ["#377eb8", "#ff7f00", "#4daf4a", "#f781bf", "#a65628", "#984ea3"]
)
bounds = np.arange(-0.5, n_clusters, 1)
norm = mcolors.BoundaryNorm(bounds, cmap.N)

# Create a figure and axis with the correct map projection
fig, ax = plt.subplots(
    figsize=(6, 4), dpi=200, subplot_kw={"projection": ccrs.PlateCarree()}
)

# Plot the K-means classification results on the map
im = ax.pcolormesh(
    longitudes,
    latitudes,
    cluster_labels,
    cmap=cmap,
    norm=norm,
    transform=ccrs.PlateCarree(),
)

# Add geographic features for context
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle=":")
ax.add_feature(cfeature.STATES, linestyle="--")

# Add gridlines
ax.gridlines(draw_labels=True)

# Set the extent to zoom in to the specified region
ax.set_extent([-95, -85, 27, 33], crs=ccrs.PlateCarree())

# Add color bar with labels
cbar = plt.colorbar(
    im,
    ax=ax,
    orientation="vertical",
    pad=0.02,
    fraction=0.05,
    ticks=np.arange(n_clusters),
)
cbar.ax.set_yticklabels([f"Class {i+1}" for i in range(n_clusters)])
cbar.set_label("Water Types", rotation=270, labelpad=20)

# Add title
ax.set_title("Water Type Classification using K-means", fontsize=16)

# Show the plot
plt.show()

# Assume 'cluster_labels' contains the K-means classification results reshaped to the original spatial dimensions # Also assume that 'da' has the original latitude and longitude data latitudes = da.coords["latitude"].values longitudes = da.coords["longitude"].values # Create a custom discrete color map for K-means clusters cmap = mcolors.ListedColormap( ["#377eb8", "#ff7f00", "#4daf4a", "#f781bf", "#a65628", "#984ea3"] ) bounds = np.arange(-0.5, n_clusters, 1) norm = mcolors.BoundaryNorm(bounds, cmap.N) # Create a figure and axis with the correct map projection fig, ax = plt.subplots( figsize=(6, 4), dpi=200, subplot_kw={"projection": ccrs.PlateCarree()} ) # Plot the K-means classification results on the map im = ax.pcolormesh( longitudes, latitudes, cluster_labels, cmap=cmap, norm=norm, transform=ccrs.PlateCarree(), ) # Add geographic features for context ax.add_feature(cfeature.COASTLINE) ax.add_feature(cfeature.BORDERS, linestyle=":") ax.add_feature(cfeature.STATES, linestyle="--") # Add gridlines ax.gridlines(draw_labels=True) # Set the extent to zoom in to the specified region ax.set_extent([-95, -85, 27, 33], crs=ccrs.PlateCarree()) # Add color bar with labels cbar = plt.colorbar( im, ax=ax, orientation="vertical", pad=0.02, fraction=0.05, ticks=np.arange(n_clusters), ) cbar.ax.set_yticklabels([f"Class {i+1}" for i in range(n_clusters)]) cbar.set_label("Water Types", rotation=270, labelpad=20) # Add title ax.set_title("Water Type Classification using K-means", fontsize=16) # Show the plot plt.show()

Keans applied to selected pixels¶

In [21]:

# Assume 'cluster_labels' contains the K-means classification results reshaped to the original spatial dimensions
# Also assume that 'da' has the original latitude and longitude data
latitudes = da.coords["latitude"].values
longitudes = da.coords["longitude"].values

# Filtering condition based on wavelength values
filter_condition = (
    (da.sel(wavelength=650) > da.sel(wavelength=620))
    & (da.sel(wavelength=701) > da.sel(wavelength=681))
    & (da.sel(wavelength=701) > da.sel(wavelength=450))
)

# Apply the filtering condition to the K-means classification results
filtered_cluster_labels = np.where(filter_condition, cluster_labels, np.nan)

# Create a custom discrete color map for K-means clusters
cmap = mcolors.ListedColormap(
    ["#e41a1c", "#377eb8", "#4daf4a", "#f781bf", "#a65628", "#984ea3"]
)
bounds = np.arange(-0.5, n_clusters, 1)
norm = mcolors.BoundaryNorm(bounds, cmap.N)

# Create a figure and axis with the correct map projection
fig, ax = plt.subplots(
    figsize=(6, 4), dpi=200, subplot_kw={"projection": ccrs.PlateCarree()}
)

# Plot the filtered K-means classification results on the map
im = ax.pcolormesh(
    longitudes,
    latitudes,
    filtered_cluster_labels,
    cmap=cmap,
    norm=norm,
    transform=ccrs.PlateCarree(),
)

# Add geographic features for context
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle=":")
ax.add_feature(cfeature.STATES, linestyle="--")

# Add gridlines
ax.gridlines(draw_labels=True)

# Set the extent to zoom in to the specified region
ax.set_extent([-95, -85, 27, 33], crs=ccrs.PlateCarree())

# Add color bar with labels
cbar = plt.colorbar(
    im,
    ax=ax,
    orientation="vertical",
    pad=0.02,
    fraction=0.05,
    ticks=np.arange(n_clusters),
)
cbar.ax.set_yticklabels([f"Class {i+1}" for i in range(n_clusters)])
cbar.set_label("Water Types", rotation=270, labelpad=20)

# Add title
ax.set_title("Water Type Classification using K-means", fontsize=16)

# Show the plot
plt.show()

# Assume 'cluster_labels' contains the K-means classification results reshaped to the original spatial dimensions # Also assume that 'da' has the original latitude and longitude data latitudes = da.coords["latitude"].values longitudes = da.coords["longitude"].values # Filtering condition based on wavelength values filter_condition = ( (da.sel(wavelength=650) > da.sel(wavelength=620)) & (da.sel(wavelength=701) > da.sel(wavelength=681)) & (da.sel(wavelength=701) > da.sel(wavelength=450)) ) # Apply the filtering condition to the K-means classification results filtered_cluster_labels = np.where(filter_condition, cluster_labels, np.nan) # Create a custom discrete color map for K-means clusters cmap = mcolors.ListedColormap( ["#e41a1c", "#377eb8", "#4daf4a", "#f781bf", "#a65628", "#984ea3"] ) bounds = np.arange(-0.5, n_clusters, 1) norm = mcolors.BoundaryNorm(bounds, cmap.N) # Create a figure and axis with the correct map projection fig, ax = plt.subplots( figsize=(6, 4), dpi=200, subplot_kw={"projection": ccrs.PlateCarree()} ) # Plot the filtered K-means classification results on the map im = ax.pcolormesh( longitudes, latitudes, filtered_cluster_labels, cmap=cmap, norm=norm, transform=ccrs.PlateCarree(), ) # Add geographic features for context ax.add_feature(cfeature.COASTLINE) ax.add_feature(cfeature.BORDERS, linestyle=":") ax.add_feature(cfeature.STATES, linestyle="--") # Add gridlines ax.gridlines(draw_labels=True) # Set the extent to zoom in to the specified region ax.set_extent([-95, -85, 27, 33], crs=ccrs.PlateCarree()) # Add color bar with labels cbar = plt.colorbar( im, ax=ax, orientation="vertical", pad=0.02, fraction=0.05, ticks=np.arange(n_clusters), ) cbar.ax.set_yticklabels([f"Class {i+1}" for i in range(n_clusters)]) cbar.set_label("Water Types", rotation=270, labelpad=20) # Add title ax.set_title("Water Type Classification using K-means", fontsize=16) # Show the plot plt.show()

PCA + SAM applied to whole image¶

In [22]:

from sklearn.decomposition import PCA

from sklearn.decomposition import PCA

In [23]:

# Check the dimensions and coordinates
print(da.dims)
print(da.coords)  # This will show the available spectral bands

# Check the dimensions and coordinates print(da.dims) print(da.coords) # This will show the available spectral bands

('latitude', 'longitude', 'wavelength')
Coordinates:
    longitude   (latitude, longitude) float32 9MB -92.88 -92.84 ... -70.54
    latitude    (latitude, longitude) float32 9MB 15.71 15.72 ... 43.91 43.92
  * wavelength  (wavelength) float64 1kB 339.0 341.0 344.0 ... 714.0 717.0 719.0

In [24]:

# Reshape data to (n_pixels, n_bands)
reshaped_data = da.values.reshape(-1, da.shape[-1])

# Handle NaNs by removing them
reshaped_data_no_nan = reshaped_data[~np.isnan(reshaped_data).any(axis=1)]

# Reshape data to (n_pixels, n_bands) reshaped_data = da.values.reshape(-1, da.shape[-1]) # Handle NaNs by removing them reshaped_data_no_nan = reshaped_data[~np.isnan(reshaped_data).any(axis=1)]

In [25]:

# Apply PCA to reduce dimensionality
pca = PCA(n_components=3)
pca_data = pca.fit_transform(reshaped_data_no_nan)

# Visualize PCA components to manually identify endmembers
plt.figure(figsize=(10, 8))
plt.scatter(pca_data[:, 0], pca_data[:, 1], c="blue", s=1)
plt.title("PCA of Spectral Data")
plt.xlabel("Principal Component 1")
plt.ylabel("Principal Component 2")
plt.show()

# Apply PCA to reduce dimensionality pca = PCA(n_components=3) pca_data = pca.fit_transform(reshaped_data_no_nan) # Visualize PCA components to manually identify endmembers plt.figure(figsize=(10, 8)) plt.scatter(pca_data[:, 0], pca_data[:, 1], c="blue", s=1) plt.title("PCA of Spectral Data") plt.xlabel("Principal Component 1") plt.ylabel("Principal Component 2") plt.show()

In [26]:

# Apply K-means to find clusters representing endmembers
n_clusters = 6  # Number of endmembers you want to find
kmeans = KMeans(n_clusters=n_clusters, random_state=0)
kmeans.fit(pca_data)

# The cluster centers in the original spectral space are your endmembers
endmembers = pca.inverse_transform(kmeans.cluster_centers_)

# Apply K-means to find clusters representing endmembers n_clusters = 6 # Number of endmembers you want to find kmeans = KMeans(n_clusters=n_clusters, random_state=0) kmeans.fit(pca_data) # The cluster centers in the original spectral space are your endmembers endmembers = pca.inverse_transform(kmeans.cluster_centers_)

In [27]:

def spectral_angle_mapper(pixel, reference):
    norm_pixel = np.linalg.norm(pixel)
    norm_reference = np.linalg.norm(reference)
    cos_theta = np.dot(pixel, reference) / (norm_pixel * norm_reference)
    angle = np.arccos(np.clip(cos_theta, -1, 1))
    return angle


# Apply SAM for each pixel and each endmember
angles = np.zeros((reshaped_data_no_nan.shape[0], endmembers.shape[0]))

for i in range(reshaped_data_no_nan.shape[0]):
    for j in range(endmembers.shape[0]):
        angles[i, j] = spectral_angle_mapper(
            reshaped_data_no_nan[i, :], endmembers[j, :]
        )

# Find the minimum angle (best match) for each pixel
best_match = np.argmin(angles, axis=1)

# Reshape best_match back to the original spatial dimensions
original_shape = da.shape[:-1]  # Get the spatial dimensions
best_match_full = np.full(reshaped_data.shape[0], np.nan)
best_match_full[~np.isnan(reshaped_data).any(axis=1)] = best_match
best_match_full = best_match_full.reshape(original_shape)

def spectral_angle_mapper(pixel, reference): norm_pixel = np.linalg.norm(pixel) norm_reference = np.linalg.norm(reference) cos_theta = np.dot(pixel, reference) / (norm_pixel * norm_reference) angle = np.arccos(np.clip(cos_theta, -1, 1)) return angle # Apply SAM for each pixel and each endmember angles = np.zeros((reshaped_data_no_nan.shape[0], endmembers.shape[0])) for i in range(reshaped_data_no_nan.shape[0]): for j in range(endmembers.shape[0]): angles[i, j] = spectral_angle_mapper( reshaped_data_no_nan[i, :], endmembers[j, :] ) # Find the minimum angle (best match) for each pixel best_match = np.argmin(angles, axis=1) # Reshape best_match back to the original spatial dimensions original_shape = da.shape[:-1] # Get the spatial dimensions best_match_full = np.full(reshaped_data.shape[0], np.nan) best_match_full[~np.isnan(reshaped_data).any(axis=1)] = best_match best_match_full = best_match_full.reshape(original_shape)

In [28]:

# Assume 'best_match_full' contains the SAM classification results reshaped to the original spatial dimensions
# Also assume that 'da' has the original latitude and longitude data
latitudes = da.coords["latitude"].values
longitudes = da.coords["longitude"].values

# Create a custom discrete color map
cmap = mcolors.ListedColormap(
    ["#377eb8", "#ff7f00", "#4daf4a", "#f781bf", "#a65628", "#984ea3"]
)
bounds = np.arange(-0.5, n_clusters, 1)
norm = mcolors.BoundaryNorm(bounds, cmap.N)

# Create a figure and axis with the correct map projection
fig, ax = plt.subplots(figsize=(10, 8), subplot_kw={"projection": ccrs.PlateCarree()})

# Plot the SAM classification results
im = ax.pcolormesh(
    longitudes,
    latitudes,
    best_match_full,
    cmap=cmap,
    norm=norm,
    transform=ccrs.PlateCarree(),
)

# Add geographic features for context
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle=":")
ax.add_feature(cfeature.STATES, linestyle="--")

# Add gridlines
ax.gridlines(draw_labels=True)

# Add color bar with labels
cbar = plt.colorbar(
    im,
    ax=ax,
    orientation="vertical",
    pad=0.02,
    fraction=0.05,
    ticks=np.arange(n_clusters),
)
cbar.ax.set_yticklabels([f"Class {i+1}" for i in range(n_clusters)])
cbar.set_label("Water Types", rotation=270, labelpad=20)

# Add title
ax.set_title("Spectral Angle Mapper (SAM) Water Type Classification", fontsize=16)

# Show the plot
plt.show()

# Assume 'best_match_full' contains the SAM classification results reshaped to the original spatial dimensions # Also assume that 'da' has the original latitude and longitude data latitudes = da.coords["latitude"].values longitudes = da.coords["longitude"].values # Create a custom discrete color map cmap = mcolors.ListedColormap( ["#377eb8", "#ff7f00", "#4daf4a", "#f781bf", "#a65628", "#984ea3"] ) bounds = np.arange(-0.5, n_clusters, 1) norm = mcolors.BoundaryNorm(bounds, cmap.N) # Create a figure and axis with the correct map projection fig, ax = plt.subplots(figsize=(10, 8), subplot_kw={"projection": ccrs.PlateCarree()}) # Plot the SAM classification results im = ax.pcolormesh( longitudes, latitudes, best_match_full, cmap=cmap, norm=norm, transform=ccrs.PlateCarree(), ) # Add geographic features for context ax.add_feature(cfeature.COASTLINE) ax.add_feature(cfeature.BORDERS, linestyle=":") ax.add_feature(cfeature.STATES, linestyle="--") # Add gridlines ax.gridlines(draw_labels=True) # Add color bar with labels cbar = plt.colorbar( im, ax=ax, orientation="vertical", pad=0.02, fraction=0.05, ticks=np.arange(n_clusters), ) cbar.ax.set_yticklabels([f"Class {i+1}" for i in range(n_clusters)]) cbar.set_label("Water Types", rotation=270, labelpad=20) # Add title ax.set_title("Spectral Angle Mapper (SAM) Water Type Classification", fontsize=16) # Show the plot plt.show()

In [29]:

# Assume 'best_match_full' contains the SAM classification results reshaped to the original spatial dimensions
# Also assume that 'da' has the original latitude and longitude data
latitudes = da.coords["latitude"].values
longitudes = da.coords["longitude"].values

# Create a custom discrete color map
cmap = mcolors.ListedColormap(
    ["#377eb8", "#ff7f00", "#4daf4a", "#f781bf", "#a65628", "#984ea3"]
)
bounds = np.arange(-0.5, n_clusters, 1)
norm = mcolors.BoundaryNorm(bounds, cmap.N)

# Create a figure and axis with the correct map projection
fig, ax = plt.subplots(figsize=(10, 8), subplot_kw={"projection": ccrs.PlateCarree()})

# Plot the SAM classification results
im = ax.pcolormesh(
    longitudes,
    latitudes,
    best_match_full,
    cmap=cmap,
    norm=norm,
    transform=ccrs.PlateCarree(),
)

# Add geographic features for context
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle=":")
ax.add_feature(cfeature.STATES, linestyle="--")

# Add gridlines
ax.gridlines(draw_labels=True)

# Set the extent to zoom in to the specified region
ax.set_extent([-95, -85, 27, 33], crs=ccrs.PlateCarree())

# Add color bar with labels
cbar = plt.colorbar(
    im,
    ax=ax,
    orientation="vertical",
    pad=0.02,
    fraction=0.05,
    ticks=np.arange(n_clusters),
)
cbar.ax.set_yticklabels([f"Class {i+1}" for i in range(n_clusters)])
cbar.set_label("Water Types", rotation=270, labelpad=20)

# Add title
ax.set_title("Spectral Angle Mapper (SAM) Water Type Classification", fontsize=16)

# Show the plot
plt.show()

# Assume 'best_match_full' contains the SAM classification results reshaped to the original spatial dimensions # Also assume that 'da' has the original latitude and longitude data latitudes = da.coords["latitude"].values longitudes = da.coords["longitude"].values # Create a custom discrete color map cmap = mcolors.ListedColormap( ["#377eb8", "#ff7f00", "#4daf4a", "#f781bf", "#a65628", "#984ea3"] ) bounds = np.arange(-0.5, n_clusters, 1) norm = mcolors.BoundaryNorm(bounds, cmap.N) # Create a figure and axis with the correct map projection fig, ax = plt.subplots(figsize=(10, 8), subplot_kw={"projection": ccrs.PlateCarree()}) # Plot the SAM classification results im = ax.pcolormesh( longitudes, latitudes, best_match_full, cmap=cmap, norm=norm, transform=ccrs.PlateCarree(), ) # Add geographic features for context ax.add_feature(cfeature.COASTLINE) ax.add_feature(cfeature.BORDERS, linestyle=":") ax.add_feature(cfeature.STATES, linestyle="--") # Add gridlines ax.gridlines(draw_labels=True) # Set the extent to zoom in to the specified region ax.set_extent([-95, -85, 27, 33], crs=ccrs.PlateCarree()) # Add color bar with labels cbar = plt.colorbar( im, ax=ax, orientation="vertical", pad=0.02, fraction=0.05, ticks=np.arange(n_clusters), ) cbar.ax.set_yticklabels([f"Class {i+1}" for i in range(n_clusters)]) cbar.set_label("Water Types", rotation=270, labelpad=20) # Add title ax.set_title("Spectral Angle Mapper (SAM) Water Type Classification", fontsize=16) # Show the plot plt.show()

In [30]:

# Assume 'best_match_full' contains the SAM classification results reshaped to the original spatial dimensions
# Also assume that 'da' has the original latitude and longitude data
latitudes = da.coords["latitude"].values
longitudes = da.coords["longitude"].values

# Filtering condition based on wavelength values
filter_condition = (
    (da.sel(wavelength=650) > da.sel(wavelength=620))
    & (da.sel(wavelength=701) > da.sel(wavelength=681))
    & (da.sel(wavelength=701) > da.sel(wavelength=450))
)

# Apply the filtering condition to the SAM classification results
filtered_best_match_full = np.where(filter_condition, best_match_full, np.nan)

# Create a custom discrete color map
cmap = mcolors.ListedColormap(
    ["#e41a1c", "#377eb8", "#4daf4a", "#f781bf", "#a65628", "#984ea3"]
)
bounds = np.arange(-0.5, n_clusters, 1)
norm = mcolors.BoundaryNorm(bounds, cmap.N)

# Create a figure and axis with the correct map projection
fig, ax = plt.subplots(figsize=(10, 8), subplot_kw={"projection": ccrs.PlateCarree()})

# Plot the filtered SAM classification results
im = ax.pcolormesh(
    longitudes,
    latitudes,
    filtered_best_match_full,
    cmap=cmap,
    norm=norm,
    transform=ccrs.PlateCarree(),
)

# Add geographic features for context
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle=":")
ax.add_feature(cfeature.STATES, linestyle="--")

# Add gridlines
ax.gridlines(draw_labels=True)

# Set the extent to zoom in to the specified region
ax.set_extent([-95, -85, 27, 33], crs=ccrs.PlateCarree())

# Add color bar with labels
cbar = plt.colorbar(
    im,
    ax=ax,
    orientation="vertical",
    pad=0.02,
    fraction=0.05,
    ticks=np.arange(n_clusters),
)
cbar.ax.set_yticklabels([f"Class {i+1}" for i in range(n_clusters)])
cbar.set_label("Water Types", rotation=270, labelpad=20)

# Add title
ax.set_title("PCA + Spectral Angle Mapper (SAM) Water Type Classification", fontsize=16)

# Show the plot
plt.show()

# Assume 'best_match_full' contains the SAM classification results reshaped to the original spatial dimensions # Also assume that 'da' has the original latitude and longitude data latitudes = da.coords["latitude"].values longitudes = da.coords["longitude"].values # Filtering condition based on wavelength values filter_condition = ( (da.sel(wavelength=650) > da.sel(wavelength=620)) & (da.sel(wavelength=701) > da.sel(wavelength=681)) & (da.sel(wavelength=701) > da.sel(wavelength=450)) ) # Apply the filtering condition to the SAM classification results filtered_best_match_full = np.where(filter_condition, best_match_full, np.nan) # Create a custom discrete color map cmap = mcolors.ListedColormap( ["#e41a1c", "#377eb8", "#4daf4a", "#f781bf", "#a65628", "#984ea3"] ) bounds = np.arange(-0.5, n_clusters, 1) norm = mcolors.BoundaryNorm(bounds, cmap.N) # Create a figure and axis with the correct map projection fig, ax = plt.subplots(figsize=(10, 8), subplot_kw={"projection": ccrs.PlateCarree()}) # Plot the filtered SAM classification results im = ax.pcolormesh( longitudes, latitudes, filtered_best_match_full, cmap=cmap, norm=norm, transform=ccrs.PlateCarree(), ) # Add geographic features for context ax.add_feature(cfeature.COASTLINE) ax.add_feature(cfeature.BORDERS, linestyle=":") ax.add_feature(cfeature.STATES, linestyle="--") # Add gridlines ax.gridlines(draw_labels=True) # Set the extent to zoom in to the specified region ax.set_extent([-95, -85, 27, 33], crs=ccrs.PlateCarree()) # Add color bar with labels cbar = plt.colorbar( im, ax=ax, orientation="vertical", pad=0.02, fraction=0.05, ticks=np.arange(n_clusters), ) cbar.ax.set_yticklabels([f"Class {i+1}" for i in range(n_clusters)]) cbar.set_label("Water Types", rotation=270, labelpad=20) # Add title ax.set_title("PCA + Spectral Angle Mapper (SAM) Water Type Classification", fontsize=16) # Show the plot plt.show()

Spectral Angle Mapper (SAM)¶

In [31]:

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
from scipy.interpolate import interp1d

import numpy as np import pandas as pd import matplotlib.pyplot as plt import matplotlib.colors as mcolors from scipy.interpolate import interp1d

In [32]:

# Load the netCDF data
file_path = "data/PACE_OCI.20240730T181157.L2.OC_AOP.V2_0.NRT.nc"
dataset = hypercoast.read_pace(file_path)
da = dataset["Rrs"]  # Assuming 'Rrs' contains the reflectance data

# Load the netCDF data file_path = "data/PACE_OCI.20240730T181157.L2.OC_AOP.V2_0.NRT.nc" dataset = hypercoast.read_pace(file_path) da = dataset["Rrs"] # Assuming 'Rrs' contains the reflectance data

In [33]:

# Extract PACE wavelengths
pace_wavelengths = da["wavelength"].values
print(pace_wavelengths)

# Extract PACE wavelengths pace_wavelengths = da["wavelength"].values print(pace_wavelengths)

[339. 341. 344. 346. 348. 351. 353. 356. 358. 361. 363. 366. 368. 371.
 373. 375. 378. 380. 383. 385. 388. 390. 393. 395. 398. 400. 403. 405.
 408. 410. 413. 415. 418. 420. 422. 425. 427. 430. 432. 435. 437. 440.
 442. 445. 447. 450. 452. 455. 457. 460. 462. 465. 467. 470. 472. 475.
 477. 480. 482. 485. 487. 490. 492. 495. 497. 500. 502. 505. 507. 510.
 512. 515. 517. 520. 522. 525. 527. 530. 532. 535. 537. 540. 542. 545.
 547. 550. 553. 555. 558. 560. 563. 565. 568. 570. 573. 575. 578. 580.
 583. 586. 588. 591. 593. 596. 598. 601. 603. 605. 608. 610. 613. 615.
 618. 620. 623. 625. 627. 630. 632. 635. 637. 640. 641. 642. 643. 645.
 646. 647. 648. 650. 651. 652. 653. 655. 656. 657. 658. 660. 661. 662.
 663. 665. 666. 667. 668. 670. 671. 672. 673. 675. 676. 677. 678. 679.
 681. 682. 683. 684. 686. 687. 688. 689. 691. 692. 693. 694. 696. 697.
 698. 699. 701. 702. 703. 704. 706. 707. 708. 709. 711. 712. 713. 714.
 717. 719.]

In [34]:

# Function to load and resample a single CSV spectral library file
def load_and_resample_spectral_library(csv_path, target_wavelengths):
    df = pd.read_csv(csv_path)
    original_wavelengths = df.iloc[:, 0].values  # First column is wavelength
    spectra_values = df.iloc[:, 1].values  # Second column is spectral values

    # Interpolation function
    interp_func = interp1d(
        original_wavelengths, spectra_values, kind="linear", fill_value="extrapolate"
    )

    # Resample to the target (PACE) wavelengths
    resampled_spectra = interp_func(target_wavelengths)

    return resampled_spectra

# Function to load and resample a single CSV spectral library file def load_and_resample_spectral_library(csv_path, target_wavelengths): df = pd.read_csv(csv_path) original_wavelengths = df.iloc[:, 0].values # First column is wavelength spectra_values = df.iloc[:, 1].values # Second column is spectral values # Interpolation function interp_func = interp1d( original_wavelengths, spectra_values, kind="linear", fill_value="extrapolate" ) # Resample to the target (PACE) wavelengths resampled_spectra = interp_func(target_wavelengths) return resampled_spectra

Download the SAM spectral library:

In [35]:

url = "https://github.com/opengeos/datasets/releases/download/hypercoast/SAM_spectral_library.zip"
hypercoast.download_file(url)

url = "https://github.com/opengeos/datasets/releases/download/hypercoast/SAM_spectral_library.zip" hypercoast.download_file(url)

Out[35]:

'/home/runner/work/HyperCoast/HyperCoast/docs/examples/SAM_spectral_library.zip'

In [36]:

# Load and resample all 6 endmembers
endmember_paths = [
    "./SAM_spectral_library/Dataset_1.csv",
    "./SAM_spectral_library/Dataset_2.csv",
    "./SAM_spectral_library/Dataset_3.csv",
    "./SAM_spectral_library/Dataset_4.csv",
    "./SAM_spectral_library/Dataset_5.csv",
    "./SAM_spectral_library/Dataset_6.csv",
]

endmembers = np.array(
    [
        load_and_resample_spectral_library(path, pace_wavelengths)
        for path in endmember_paths
    ]
)

# Load and resample all 6 endmembers endmember_paths = [ "./SAM_spectral_library/Dataset_1.csv", "./SAM_spectral_library/Dataset_2.csv", "./SAM_spectral_library/Dataset_3.csv", "./SAM_spectral_library/Dataset_4.csv", "./SAM_spectral_library/Dataset_5.csv", "./SAM_spectral_library/Dataset_6.csv", ] endmembers = np.array( [ load_and_resample_spectral_library(path, pace_wavelengths) for path in endmember_paths ] )

In [37]:

print(endmembers)

print(endmembers)

[[0.00871842 0.00865469 0.00855908 ... 0.00442423 0.00389996 0.00375134]
 [0.01088283 0.01082291 0.01073303 ... 0.01052424 0.00998495 0.00957979]
 [0.0169974  0.0168173  0.01654716 ... 0.01991858 0.01881549 0.01798281]
 [0.0227546  0.022496   0.0221081  ... 0.0312176  0.03006733 0.02898674]
 [0.02859302 0.02826977 0.02778489 ... 0.05700807 0.05723728 0.05712533]
 [0.03042756 0.03000733 0.02937699 ... 0.07612622 0.08101502 0.08450761]]

In [38]:

# Plot sample spectra from the CSV files and their resampled versions
def plot_sample_spectra(csv_paths, pace_wavelengths):
    plt.figure(figsize=(14, 8))

    for i, csv_path in enumerate(csv_paths):
        df = pd.read_csv(csv_path)
        original_wavelengths = df.iloc[:, 0].values
        spectra_values = df.iloc[:, 1].values
        resampled_spectra = load_and_resample_spectral_library(
            csv_path, pace_wavelengths
        )

        plt.plot(
            original_wavelengths,
            spectra_values,
            label=f"Original Spectra {i+1}",
            linestyle="--",
        )
        plt.plot(pace_wavelengths, resampled_spectra, label=f"Resampled Spectra {i+1}")

    plt.xlabel("Wavelength (nm)")
    plt.ylabel("Spectral Reflectance")
    plt.title("Comparison of Original and Resampled Spectra")
    plt.legend()
    plt.grid(True)
    plt.show()


plot_sample_spectra(endmember_paths, pace_wavelengths)

# Plot sample spectra from the CSV files and their resampled versions def plot_sample_spectra(csv_paths, pace_wavelengths): plt.figure(figsize=(14, 8)) for i, csv_path in enumerate(csv_paths): df = pd.read_csv(csv_path) original_wavelengths = df.iloc[:, 0].values spectra_values = df.iloc[:, 1].values resampled_spectra = load_and_resample_spectral_library( csv_path, pace_wavelengths ) plt.plot( original_wavelengths, spectra_values, label=f"Original Spectra {i+1}", linestyle="--", ) plt.plot(pace_wavelengths, resampled_spectra, label=f"Resampled Spectra {i+1}") plt.xlabel("Wavelength (nm)") plt.ylabel("Spectral Reflectance") plt.title("Comparison of Original and Resampled Spectra") plt.legend() plt.grid(True) plt.show() plot_sample_spectra(endmember_paths, pace_wavelengths)

In [39]:

# Function to calculate spectral angle
def spectral_angle_mapper(pixel, reference):
    norm_pixel = np.linalg.norm(pixel)
    norm_reference = np.linalg.norm(reference)
    cos_theta = np.dot(pixel, reference) / (norm_pixel * norm_reference)
    angle = np.arccos(np.clip(cos_theta, -1, 1))
    return angle

# Function to calculate spectral angle def spectral_angle_mapper(pixel, reference): norm_pixel = np.linalg.norm(pixel) norm_reference = np.linalg.norm(reference) cos_theta = np.dot(pixel, reference) / (norm_pixel * norm_reference) angle = np.arccos(np.clip(cos_theta, -1, 1)) return angle

In [40]:

# Reshape data to (n_pixels, n_bands)
reshaped_data = da.values.reshape(-1, da.shape[-1])

# Reshape data to (n_pixels, n_bands) reshaped_data = da.values.reshape(-1, da.shape[-1])

In [41]:

# Apply SAM for each pixel and each endmember
angles = np.zeros((reshaped_data.shape[0], endmembers.shape[0]))

for i in range(reshaped_data.shape[0]):
    for j in range(endmembers.shape[0]):
        angles[i, j] = spectral_angle_mapper(reshaped_data[i, :], endmembers[j, :])

# Apply SAM for each pixel and each endmember angles = np.zeros((reshaped_data.shape[0], endmembers.shape[0])) for i in range(reshaped_data.shape[0]): for j in range(endmembers.shape[0]): angles[i, j] = spectral_angle_mapper(reshaped_data[i, :], endmembers[j, :])

In [42]:

# Find the minimum angle (best match) for each pixel
best_match = np.argmin(angles, axis=1)

# Find the minimum angle (best match) for each pixel best_match = np.argmin(angles, axis=1)

In [43]:

# Reshape best_match back to the original spatial dimensions
best_match = best_match.reshape(da.shape[:-1])

# Reshape best_match back to the original spatial dimensions best_match = best_match.reshape(da.shape[:-1])

In [44]:

# Assume 'best_match' contains the SAM classification results reshaped to the original spatial dimensions
# Also assume that 'da' has the original latitude and longitude data
latitudes = da.coords["latitude"].values
longitudes = da.coords["longitude"].values

# Create a custom discrete color map
cmap = mcolors.ListedColormap(
    ["#377eb8", "#e41a1c", "#4daf4a", "#f781bf", "#a65628", "#984ea3"]
)
bounds = np.arange(-0.5, len(endmembers), 1)
norm = mcolors.BoundaryNorm(bounds, cmap.N)

# Create a figure and axis with the correct map projection
fig, ax = plt.subplots(figsize=(12, 10), subplot_kw={"projection": ccrs.PlateCarree()})

# Plot the SAM classification results
im = ax.pcolormesh(
    longitudes,
    latitudes,
    best_match,
    cmap=cmap,
    norm=norm,
    transform=ccrs.PlateCarree(),
)

# Add geographic features for context
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle=":")
ax.add_feature(cfeature.STATES, linestyle="--")

# Adding axis labels
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")

# Adding a title
ax.set_title("Spectral Angle Mapper (SAM) Water Type Classification", fontsize=16)

# Adding a color bar with discrete values
cbar = plt.colorbar(
    im,
    ax=ax,
    orientation="vertical",
    pad=0.02,
    fraction=0.05,
    ticks=np.arange(len(endmembers)),
)
cbar.ax.set_yticklabels([f"Class {i+1}" for i in range(len(endmembers))])
cbar.set_label("Water Types", rotation=270, labelpad=20)

# Adding gridlines
ax.gridlines(draw_labels=True, linestyle="--", linewidth=0.5)

# Set the extent to zoom in to the specified region (adjust as needed)
ax.set_extent([-95, -85, 27, 33], crs=ccrs.PlateCarree())

# Show the plot
plt.show()

# Assume 'best_match' contains the SAM classification results reshaped to the original spatial dimensions # Also assume that 'da' has the original latitude and longitude data latitudes = da.coords["latitude"].values longitudes = da.coords["longitude"].values # Create a custom discrete color map cmap = mcolors.ListedColormap( ["#377eb8", "#e41a1c", "#4daf4a", "#f781bf", "#a65628", "#984ea3"] ) bounds = np.arange(-0.5, len(endmembers), 1) norm = mcolors.BoundaryNorm(bounds, cmap.N) # Create a figure and axis with the correct map projection fig, ax = plt.subplots(figsize=(12, 10), subplot_kw={"projection": ccrs.PlateCarree()}) # Plot the SAM classification results im = ax.pcolormesh( longitudes, latitudes, best_match, cmap=cmap, norm=norm, transform=ccrs.PlateCarree(), ) # Add geographic features for context ax.add_feature(cfeature.COASTLINE) ax.add_feature(cfeature.BORDERS, linestyle=":") ax.add_feature(cfeature.STATES, linestyle="--") # Adding axis labels ax.set_xlabel("Longitude") ax.set_ylabel("Latitude") # Adding a title ax.set_title("Spectral Angle Mapper (SAM) Water Type Classification", fontsize=16) # Adding a color bar with discrete values cbar = plt.colorbar( im, ax=ax, orientation="vertical", pad=0.02, fraction=0.05, ticks=np.arange(len(endmembers)), ) cbar.ax.set_yticklabels([f"Class {i+1}" for i in range(len(endmembers))]) cbar.set_label("Water Types", rotation=270, labelpad=20) # Adding gridlines ax.gridlines(draw_labels=True, linestyle="--", linewidth=0.5) # Set the extent to zoom in to the specified region (adjust as needed) ax.set_extent([-95, -85, 27, 33], crs=ccrs.PlateCarree()) # Show the plot plt.show()

In [45]:

# Assume 'best_match' contains the SAM classification results reshaped to the original spatial dimensions
# Also assume that 'da' has the original latitude and longitude data
latitudes = da.coords["latitude"].values
longitudes = da.coords["longitude"].values

# Extract specific wavelengths for the conditions
wavelength_450 = da.sel(wavelength=450).values
wavelength_620 = da.sel(wavelength=620).values
wavelength_650 = da.sel(wavelength=650).values
wavelength_681 = da.sel(wavelength=681).values
wavelength_701 = da.sel(wavelength=701).values

# Apply the condition to filter the pixels
condition = (
    (wavelength_650 > wavelength_620)
    & (wavelength_701 > wavelength_681)
    & (wavelength_701 > wavelength_450)
)

# Filter the best_match data based on the condition
filtered_best_match = np.where(condition, best_match, np.nan)

# Create a custom discrete color map
cmap = mcolors.ListedColormap(
    ["#377eb8", "#e41a1c", "#4daf4a", "#f781bf", "#a65628", "#984ea3"]
)
bounds = np.arange(-0.5, len(endmembers), 1)
norm = mcolors.BoundaryNorm(bounds, cmap.N)

# Create a figure and axis with the correct map projection
fig, ax = plt.subplots(figsize=(12, 10), subplot_kw={"projection": ccrs.PlateCarree()})

# Plot the SAM classification results for the filtered pixels
im = ax.pcolormesh(
    longitudes,
    latitudes,
    filtered_best_match,
    cmap=cmap,
    norm=norm,
    transform=ccrs.PlateCarree(),
)

# Add geographic features for context
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle=":")
ax.add_feature(cfeature.STATES, linestyle="--")

# Adding axis labels
ax.set_xlabel("Longitude")
ax.set_ylabel("Latitude")

# Adding a title
ax.set_title("SAM Water Type Classification", fontsize=16)

# Adding a color bar with discrete values
cbar = plt.colorbar(
    im,
    ax=ax,
    orientation="vertical",
    pad=0.02,
    fraction=0.05,
    ticks=np.arange(len(endmembers)),
)
cbar.ax.set_yticklabels([f"Class {i+1}" for i in range(len(endmembers))])
cbar.set_label("Water Types", rotation=270, labelpad=20)

# Adding gridlines
ax.gridlines(draw_labels=True, linestyle="--", linewidth=0.5)

# Set the extent to zoom in to the specified region (adjust as needed)
ax.set_extent([-95, -85, 27, 33], crs=ccrs.PlateCarree())

# Show the plot
plt.show()

# Assume 'best_match' contains the SAM classification results reshaped to the original spatial dimensions # Also assume that 'da' has the original latitude and longitude data latitudes = da.coords["latitude"].values longitudes = da.coords["longitude"].values # Extract specific wavelengths for the conditions wavelength_450 = da.sel(wavelength=450).values wavelength_620 = da.sel(wavelength=620).values wavelength_650 = da.sel(wavelength=650).values wavelength_681 = da.sel(wavelength=681).values wavelength_701 = da.sel(wavelength=701).values # Apply the condition to filter the pixels condition = ( (wavelength_650 > wavelength_620) & (wavelength_701 > wavelength_681) & (wavelength_701 > wavelength_450) ) # Filter the best_match data based on the condition filtered_best_match = np.where(condition, best_match, np.nan) # Create a custom discrete color map cmap = mcolors.ListedColormap( ["#377eb8", "#e41a1c", "#4daf4a", "#f781bf", "#a65628", "#984ea3"] ) bounds = np.arange(-0.5, len(endmembers), 1) norm = mcolors.BoundaryNorm(bounds, cmap.N) # Create a figure and axis with the correct map projection fig, ax = plt.subplots(figsize=(12, 10), subplot_kw={"projection": ccrs.PlateCarree()}) # Plot the SAM classification results for the filtered pixels im = ax.pcolormesh( longitudes, latitudes, filtered_best_match, cmap=cmap, norm=norm, transform=ccrs.PlateCarree(), ) # Add geographic features for context ax.add_feature(cfeature.COASTLINE) ax.add_feature(cfeature.BORDERS, linestyle=":") ax.add_feature(cfeature.STATES, linestyle="--") # Adding axis labels ax.set_xlabel("Longitude") ax.set_ylabel("Latitude") # Adding a title ax.set_title("SAM Water Type Classification", fontsize=16) # Adding a color bar with discrete values cbar = plt.colorbar( im, ax=ax, orientation="vertical", pad=0.02, fraction=0.05, ticks=np.arange(len(endmembers)), ) cbar.ax.set_yticklabels([f"Class {i+1}" for i in range(len(endmembers))]) cbar.set_label("Water Types", rotation=270, labelpad=20) # Adding gridlines ax.gridlines(draw_labels=True, linestyle="--", linewidth=0.5) # Set the extent to zoom in to the specified region (adjust as needed) ax.set_extent([-95, -85, 27, 33], crs=ccrs.PlateCarree()) # Show the plot plt.show()