diff --git a/docs/examples/pace_cyano.ipynb b/docs/examples/pace_cyano.ipynb
index 24658b1..31b6161 100644
--- a/docs/examples/pace_cyano.ipynb
+++ b/docs/examples/pace_cyano.ipynb
@@ -49,7 +49,14 @@
    "outputs": [],
    "source": [
     "import earthaccess\n",
-    "import hypercoast"
+    "import hypercoast\n",
+    "from hypercoast.pace import (\n",
+    "    cyano_band_ratios,\n",
+    "    apply_kmeans,\n",
+    "    apply_pca,\n",
+    "    apply_sam,\n",
+    "    apply_sam_spectral,\n",
+    ")"
    ]
   },
   {
@@ -155,13 +162,7 @@
    },
    "outputs": [],
    "source": [
-    "da = dataset[\"Rrs\"]\n",
-    "data = (\n",
-    "    (da.sel(wavelength=650) > da.sel(wavelength=620))\n",
-    "    & (da.sel(wavelength=701) > da.sel(wavelength=681))\n",
-    "    & (da.sel(wavelength=701) > da.sel(wavelength=450))\n",
-    ")\n",
-    "# data"
+    "da = cyano_band_ratios(dataset, plot=True)"
    ]
   },
   {
@@ -175,65 +176,6 @@
     "![](https://i.imgur.com/pQP50bz.png)"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "metadata": {
-    "id": "EuKeNoNp-_JJ"
-   },
-   "source": [
-    "## Visualize the selected region based on band ratios"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "colab": {
-     "base_uri": "https://localhost:8080/",
-     "height": 540
-    },
-    "id": "be6r9fkf-_JJ",
-    "outputId": "a1964c49-6587-475a-87bd-a95db40bb35e"
-   },
-   "outputs": [],
-   "source": [
-    "import xarray as xr\n",
-    "import matplotlib.pyplot as plt\n",
-    "import cartopy.crs as ccrs\n",
-    "import cartopy.feature as cfeature\n",
-    "\n",
-    "# Create a plot\n",
-    "fig, ax = plt.subplots(figsize=(12, 6), subplot_kw={\"projection\": ccrs.PlateCarree()})\n",
-    "\n",
-    "# ax.set_extent([-93, -87, 28, 32], crs=ccrs.PlateCarree())\n",
-    "\n",
-    "# Plot the data\n",
-    "data.plot(\n",
-    "    ax=ax,\n",
-    "    transform=ccrs.PlateCarree(),\n",
-    "    cmap=\"coolwarm\",\n",
-    "    cbar_kwargs={\"label\": \"Cyano\"},\n",
-    ")\n",
-    "\n",
-    "# Add coastlines\n",
-    "ax.coastlines()\n",
-    "\n",
-    "# Add state boundaries\n",
-    "states_provinces = cfeature.NaturalEarthFeature(\n",
-    "    category=\"cultural\",\n",
-    "    name=\"admin_1_states_provinces_lines\",\n",
-    "    scale=\"50m\",\n",
-    "    facecolor=\"none\",\n",
-    ")\n",
-    "\n",
-    "ax.add_feature(states_provinces, edgecolor=\"gray\")\n",
-    "\n",
-    "# Optionally, add gridlines, labels, etc.\n",
-    "ax.gridlines(draw_labels=True)\n",
-    "\n",
-    "plt.show()"
-   ]
-  },
   {
    "cell_type": "markdown",
    "metadata": {
@@ -257,230 +199,20 @@
     "## K-means applied to the whole image"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "id": "jDrdL13RcLFm"
-   },
-   "outputs": [],
-   "source": [
-    "import numpy as np\n",
-    "from sklearn.cluster import KMeans"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "colab": {
-     "base_uri": "https://localhost:8080/"
-    },
-    "id": "8SYB8MKyZc0-",
-    "outputId": "02c93296-323b-41e9-e175-eaa18e0d19eb"
-   },
-   "outputs": [],
-   "source": [
-    "# Get the shape of the DataArray\n",
-    "print(\"Shape of da:\", da.shape)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "colab": {
-     "base_uri": "https://localhost:8080/"
-    },
-    "id": "eCH5c_uBZc28",
-    "outputId": "1301bab1-74d6-4507-ce27-b1f0ae81d943"
-   },
-   "outputs": [],
-   "source": [
-    "# Get the dimension names of the DataArray\n",
-    "print(\"Dimension names:\", da.dims)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "colab": {
-     "base_uri": "https://localhost:8080/"
-    },
-    "id": "0nzjK69bZc6d",
-    "outputId": "fd909e7d-dcc2-451c-f8aa-d67545ee7673"
-   },
-   "outputs": [],
-   "source": [
-    "# Get the size of each dimension\n",
-    "print(\"Size of each dimension:\", da.sizes)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "id": "_O_obmMXZc9Z"
-   },
-   "outputs": [],
-   "source": [
-    "reshaped_data = da.values.reshape(-1, da.shape[-1])"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "id": "1V22cNgbZdAz"
-   },
-   "outputs": [],
-   "source": [
-    "reshaped_data_no_nan = reshaped_data[~np.isnan(reshaped_data).any(axis=1)]"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "colab": {
-     "base_uri": "https://localhost:8080/",
-     "height": 75
-    },
-    "id": "70hwmWwZelJn",
-    "outputId": "b0daf9fc-3b4a-400e-e512-c40de1ae885a"
-   },
-   "outputs": [],
-   "source": [
-    "# Apply K-means clustering to classify into 5-6 water types.\n",
-    "n_clusters = 6\n",
-    "kmeans = KMeans(n_clusters=n_clusters, random_state=0)\n",
-    "kmeans.fit(reshaped_data_no_nan)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "id": "yF1Y40bhelMn"
-   },
-   "outputs": [],
-   "source": [
-    "# Initialize an array for cluster labels with NaN\n",
-    "labels = np.full(reshaped_data.shape[0], np.nan)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "id": "nIrI-WyGelPj"
-   },
-   "outputs": [],
-   "source": [
-    "# Assign the computed cluster labels to the non-NaN positions\n",
-    "labels[~np.isnan(reshaped_data).any(axis=1)] = kmeans.labels_"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "id": "6u86WCEnfnUU"
-   },
-   "outputs": [],
-   "source": [
-    "# Reshape the labels back to the original spatial dimensions\n",
-    "cluster_labels = labels.reshape(da.shape[:-1])"
-   ]
-  },
   {
    "cell_type": "code",
    "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
-    "import matplotlib.pyplot as plt\n",
-    "import matplotlib.colors as mcolors\n",
-    "import cartopy.crs as ccrs\n",
-    "import cartopy.feature as cfeature"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "colab": {
-     "base_uri": "https://localhost:8080/",
-     "height": 704
-    },
-    "id": "K7mEYr4iny5f",
-    "outputId": "394a976d-d9d7-496c-9dec-0764dc4f1bf5"
-   },
-   "outputs": [],
-   "source": [
-    "# Assume 'cluster_labels' contains the K-means classification results reshaped to the original spatial dimensions\n",
-    "# Also assume that 'da' has the original latitude and longitude data\n",
-    "latitudes = da.coords[\"latitude\"].values\n",
-    "longitudes = da.coords[\"longitude\"].values\n",
-    "\n",
-    "# Create a custom discrete color map for K-means clusters\n",
-    "cmap = mcolors.ListedColormap(\n",
-    "    [\"#377eb8\", \"#ff7f00\", \"#4daf4a\", \"#f781bf\", \"#a65628\", \"#984ea3\"]\n",
-    ")\n",
-    "bounds = np.arange(-0.5, n_clusters, 1)\n",
-    "norm = mcolors.BoundaryNorm(bounds, cmap.N)\n",
-    "\n",
-    "# Create a figure and axis with the correct map projection\n",
-    "fig, ax = plt.subplots(\n",
-    "    figsize=(6, 4), dpi=200, subplot_kw={\"projection\": ccrs.PlateCarree()}\n",
-    ")\n",
-    "\n",
-    "# Plot the K-means classification results on the map\n",
-    "im = ax.pcolormesh(\n",
-    "    longitudes,\n",
-    "    latitudes,\n",
-    "    cluster_labels,\n",
-    "    cmap=cmap,\n",
-    "    norm=norm,\n",
-    "    transform=ccrs.PlateCarree(),\n",
-    ")\n",
-    "\n",
-    "# Add geographic features for context\n",
-    "ax.add_feature(cfeature.COASTLINE)\n",
-    "ax.add_feature(cfeature.BORDERS, linestyle=\":\")\n",
-    "ax.add_feature(cfeature.STATES, linestyle=\"--\")\n",
-    "\n",
-    "# Add gridlines\n",
-    "ax.gridlines(draw_labels=True)\n",
-    "\n",
-    "# Set the extent to zoom in to the specified region\n",
-    "ax.set_extent([-95, -85, 27, 33], crs=ccrs.PlateCarree())\n",
-    "\n",
-    "# Add color bar with labels\n",
-    "cbar = plt.colorbar(\n",
-    "    im,\n",
-    "    ax=ax,\n",
-    "    orientation=\"vertical\",\n",
-    "    pad=0.02,\n",
-    "    fraction=0.05,\n",
-    "    ticks=np.arange(n_clusters),\n",
-    ")\n",
-    "cbar.ax.set_yticklabels([f\"Class {i+1}\" for i in range(n_clusters)])\n",
-    "cbar.set_label(\"Water Types\", rotation=270, labelpad=20)\n",
-    "\n",
-    "# Add title\n",
-    "ax.set_title(\"Water Type Classification using K-means\", fontsize=16)\n",
-    "\n",
-    "# Show the plot\n",
-    "plt.show()"
+    "cluster_labels, latitudes, longitudes = apply_kmeans(dataset, n_clusters=6)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Keans applied to selected pixels"
+    "## K-means applied to selected pixels"
    ]
   },
   {
@@ -491,423 +223,26 @@
    },
    "outputs": [],
    "source": [
-    "# Assume 'cluster_labels' contains the K-means classification results reshaped to the original spatial dimensions\n",
-    "# Also assume that 'da' has the original latitude and longitude data\n",
-    "latitudes = da.coords[\"latitude\"].values\n",
-    "longitudes = da.coords[\"longitude\"].values\n",
-    "\n",
-    "# Filtering condition based on wavelength values\n",
-    "filter_condition = (\n",
-    "    (da.sel(wavelength=650) > da.sel(wavelength=620))\n",
-    "    & (da.sel(wavelength=701) > da.sel(wavelength=681))\n",
-    "    & (da.sel(wavelength=701) > da.sel(wavelength=450))\n",
-    ")\n",
-    "\n",
-    "# Apply the filtering condition to the K-means classification results\n",
-    "filtered_cluster_labels = np.where(filter_condition, cluster_labels, np.nan)\n",
-    "\n",
-    "# Create a custom discrete color map for K-means clusters\n",
-    "cmap = mcolors.ListedColormap(\n",
-    "    [\"#e41a1c\", \"#377eb8\", \"#4daf4a\", \"#f781bf\", \"#a65628\", \"#984ea3\"]\n",
-    ")\n",
-    "bounds = np.arange(-0.5, n_clusters, 1)\n",
-    "norm = mcolors.BoundaryNorm(bounds, cmap.N)\n",
-    "\n",
-    "# Create a figure and axis with the correct map projection\n",
-    "fig, ax = plt.subplots(\n",
-    "    figsize=(6, 4), dpi=200, subplot_kw={\"projection\": ccrs.PlateCarree()}\n",
-    ")\n",
-    "\n",
-    "# Plot the filtered K-means classification results on the map\n",
-    "im = ax.pcolormesh(\n",
-    "    longitudes,\n",
-    "    latitudes,\n",
-    "    filtered_cluster_labels,\n",
-    "    cmap=cmap,\n",
-    "    norm=norm,\n",
-    "    transform=ccrs.PlateCarree(),\n",
-    ")\n",
-    "\n",
-    "# Add geographic features for context\n",
-    "ax.add_feature(cfeature.COASTLINE)\n",
-    "ax.add_feature(cfeature.BORDERS, linestyle=\":\")\n",
-    "ax.add_feature(cfeature.STATES, linestyle=\"--\")\n",
-    "\n",
-    "# Add gridlines\n",
-    "ax.gridlines(draw_labels=True)\n",
-    "\n",
-    "# Set the extent to zoom in to the specified region\n",
-    "ax.set_extent([-95, -85, 27, 33], crs=ccrs.PlateCarree())\n",
-    "\n",
-    "# Add color bar with labels\n",
-    "cbar = plt.colorbar(\n",
-    "    im,\n",
-    "    ax=ax,\n",
-    "    orientation=\"vertical\",\n",
-    "    pad=0.02,\n",
-    "    fraction=0.05,\n",
-    "    ticks=np.arange(n_clusters),\n",
-    ")\n",
-    "cbar.ax.set_yticklabels([f\"Class {i+1}\" for i in range(n_clusters)])\n",
-    "cbar.set_label(\"Water Types\", rotation=270, labelpad=20)\n",
-    "\n",
-    "# Add title\n",
-    "ax.set_title(\"Water Type Classification using K-means\", fontsize=16)\n",
-    "\n",
-    "# Show the plot\n",
-    "plt.show()"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "## PCA + SAM applied to whole image"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "from sklearn.decomposition import PCA"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "colab": {
-     "base_uri": "https://localhost:8080/"
-    },
-    "id": "svg0vrnxh6vp",
-    "outputId": "f7d2d17f-5035-4349-cd2d-32801f7feb17"
-   },
-   "outputs": [],
-   "source": [
-    "# Check the dimensions and coordinates\n",
-    "print(da.dims)\n",
-    "print(da.coords)  # This will show the available spectral bands"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "id": "dksBIw9KjUfS"
-   },
-   "outputs": [],
-   "source": [
-    "# Reshape data to (n_pixels, n_bands)\n",
-    "reshaped_data = da.values.reshape(-1, da.shape[-1])\n",
-    "\n",
-    "# Handle NaNs by removing them\n",
-    "reshaped_data_no_nan = reshaped_data[~np.isnan(reshaped_data).any(axis=1)]"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "colab": {
-     "base_uri": "https://localhost:8080/",
-     "height": 718
-    },
-    "id": "f_R6a2nDjUjR",
-    "outputId": "2b8e4151-e9d3-492c-edcf-1db12e505322"
-   },
-   "outputs": [],
-   "source": [
-    "# Apply PCA to reduce dimensionality\n",
-    "pca = PCA(n_components=3)\n",
-    "pca_data = pca.fit_transform(reshaped_data_no_nan)\n",
-    "\n",
-    "# Visualize PCA components to manually identify endmembers\n",
-    "plt.figure(figsize=(10, 8))\n",
-    "plt.scatter(pca_data[:, 0], pca_data[:, 1], c=\"blue\", s=1)\n",
-    "plt.title(\"PCA of Spectral Data\")\n",
-    "plt.xlabel(\"Principal Component 1\")\n",
-    "plt.ylabel(\"Principal Component 2\")\n",
-    "plt.show()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "id": "sLCiBgBajpGC"
-   },
-   "outputs": [],
-   "source": [
-    "# Apply K-means to find clusters representing endmembers\n",
-    "n_clusters = 6  # Number of endmembers you want to find\n",
-    "kmeans = KMeans(n_clusters=n_clusters, random_state=0)\n",
-    "kmeans.fit(pca_data)\n",
-    "\n",
-    "# The cluster centers in the original spectral space are your endmembers\n",
-    "endmembers = pca.inverse_transform(kmeans.cluster_centers_)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "id": "iPkI2yuRjpR9"
-   },
-   "outputs": [],
-   "source": [
-    "def spectral_angle_mapper(pixel, reference):\n",
-    "    norm_pixel = np.linalg.norm(pixel)\n",
-    "    norm_reference = np.linalg.norm(reference)\n",
-    "    cos_theta = np.dot(pixel, reference) / (norm_pixel * norm_reference)\n",
-    "    angle = np.arccos(np.clip(cos_theta, -1, 1))\n",
-    "    return angle\n",
-    "\n",
-    "\n",
-    "# Apply SAM for each pixel and each endmember\n",
-    "angles = np.zeros((reshaped_data_no_nan.shape[0], endmembers.shape[0]))\n",
-    "\n",
-    "for i in range(reshaped_data_no_nan.shape[0]):\n",
-    "    for j in range(endmembers.shape[0]):\n",
-    "        angles[i, j] = spectral_angle_mapper(\n",
-    "            reshaped_data_no_nan[i, :], endmembers[j, :]\n",
-    "        )\n",
-    "\n",
-    "# Find the minimum angle (best match) for each pixel\n",
-    "best_match = np.argmin(angles, axis=1)\n",
-    "\n",
-    "# Reshape best_match back to the original spatial dimensions\n",
-    "original_shape = da.shape[:-1]  # Get the spatial dimensions\n",
-    "best_match_full = np.full(reshaped_data.shape[0], np.nan)\n",
-    "best_match_full[~np.isnan(reshaped_data).any(axis=1)] = best_match\n",
-    "best_match_full = best_match_full.reshape(original_shape)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "colab": {
-     "base_uri": "https://localhost:8080/",
-     "height": 704
-    },
-    "id": "3EuvKbPHjUq7",
-    "outputId": "41372b92-48e2-4c47-e732-5e597e4086a0"
-   },
-   "outputs": [],
-   "source": [
-    "# Assume 'best_match_full' contains the SAM classification results reshaped to the original spatial dimensions\n",
-    "# Also assume that 'da' has the original latitude and longitude data\n",
-    "latitudes = da.coords[\"latitude\"].values\n",
-    "longitudes = da.coords[\"longitude\"].values\n",
-    "\n",
-    "# Create a custom discrete color map\n",
-    "cmap = mcolors.ListedColormap(\n",
-    "    [\"#377eb8\", \"#ff7f00\", \"#4daf4a\", \"#f781bf\", \"#a65628\", \"#984ea3\"]\n",
-    ")\n",
-    "bounds = np.arange(-0.5, n_clusters, 1)\n",
-    "norm = mcolors.BoundaryNorm(bounds, cmap.N)\n",
-    "\n",
-    "# Create a figure and axis with the correct map projection\n",
-    "fig, ax = plt.subplots(figsize=(10, 8), subplot_kw={\"projection\": ccrs.PlateCarree()})\n",
-    "\n",
-    "# Plot the SAM classification results\n",
-    "im = ax.pcolormesh(\n",
-    "    longitudes,\n",
-    "    latitudes,\n",
-    "    best_match_full,\n",
-    "    cmap=cmap,\n",
-    "    norm=norm,\n",
-    "    transform=ccrs.PlateCarree(),\n",
-    ")\n",
-    "\n",
-    "# Add geographic features for context\n",
-    "ax.add_feature(cfeature.COASTLINE)\n",
-    "ax.add_feature(cfeature.BORDERS, linestyle=\":\")\n",
-    "ax.add_feature(cfeature.STATES, linestyle=\"--\")\n",
-    "\n",
-    "# Add gridlines\n",
-    "ax.gridlines(draw_labels=True)\n",
-    "\n",
-    "# Add color bar with labels\n",
-    "cbar = plt.colorbar(\n",
-    "    im,\n",
-    "    ax=ax,\n",
-    "    orientation=\"vertical\",\n",
-    "    pad=0.02,\n",
-    "    fraction=0.05,\n",
-    "    ticks=np.arange(n_clusters),\n",
-    ")\n",
-    "cbar.ax.set_yticklabels([f\"Class {i+1}\" for i in range(n_clusters)])\n",
-    "cbar.set_label(\"Water Types\", rotation=270, labelpad=20)\n",
-    "\n",
-    "# Add title\n",
-    "ax.set_title(\"Spectral Angle Mapper (SAM) Water Type Classification\", fontsize=16)\n",
-    "\n",
-    "# Show the plot\n",
-    "plt.show()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "colab": {
-     "base_uri": "https://localhost:8080/",
-     "height": 653
-    },
-    "id": "jnJ48tcxjUve",
-    "outputId": "70585cfb-d9e4-49e7-80a7-8393365a3094"
-   },
-   "outputs": [],
-   "source": [
-    "# Assume 'best_match_full' contains the SAM classification results reshaped to the original spatial dimensions\n",
-    "# Also assume that 'da' has the original latitude and longitude data\n",
-    "latitudes = da.coords[\"latitude\"].values\n",
-    "longitudes = da.coords[\"longitude\"].values\n",
-    "\n",
-    "# Create a custom discrete color map\n",
-    "cmap = mcolors.ListedColormap(\n",
-    "    [\"#377eb8\", \"#ff7f00\", \"#4daf4a\", \"#f781bf\", \"#a65628\", \"#984ea3\"]\n",
-    ")\n",
-    "bounds = np.arange(-0.5, n_clusters, 1)\n",
-    "norm = mcolors.BoundaryNorm(bounds, cmap.N)\n",
-    "\n",
-    "# Create a figure and axis with the correct map projection\n",
-    "fig, ax = plt.subplots(figsize=(10, 8), subplot_kw={\"projection\": ccrs.PlateCarree()})\n",
-    "\n",
-    "# Plot the SAM classification results\n",
-    "im = ax.pcolormesh(\n",
-    "    longitudes,\n",
-    "    latitudes,\n",
-    "    best_match_full,\n",
-    "    cmap=cmap,\n",
-    "    norm=norm,\n",
-    "    transform=ccrs.PlateCarree(),\n",
-    ")\n",
-    "\n",
-    "# Add geographic features for context\n",
-    "ax.add_feature(cfeature.COASTLINE)\n",
-    "ax.add_feature(cfeature.BORDERS, linestyle=\":\")\n",
-    "ax.add_feature(cfeature.STATES, linestyle=\"--\")\n",
-    "\n",
-    "# Add gridlines\n",
-    "ax.gridlines(draw_labels=True)\n",
-    "\n",
-    "# Set the extent to zoom in to the specified region\n",
-    "ax.set_extent([-95, -85, 27, 33], crs=ccrs.PlateCarree())\n",
-    "\n",
-    "# Add color bar with labels\n",
-    "cbar = plt.colorbar(\n",
-    "    im,\n",
-    "    ax=ax,\n",
-    "    orientation=\"vertical\",\n",
-    "    pad=0.02,\n",
-    "    fraction=0.05,\n",
-    "    ticks=np.arange(n_clusters),\n",
-    ")\n",
-    "cbar.ax.set_yticklabels([f\"Class {i+1}\" for i in range(n_clusters)])\n",
-    "cbar.set_label(\"Water Types\", rotation=270, labelpad=20)\n",
-    "\n",
-    "# Add title\n",
-    "ax.set_title(\"Spectral Angle Mapper (SAM) Water Type Classification\", fontsize=16)\n",
-    "\n",
-    "# Show the plot\n",
-    "plt.show()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {
-    "id": "HflML7NnjUym"
-   },
-   "outputs": [],
-   "source": [
-    "# Assume 'best_match_full' contains the SAM classification results reshaped to the original spatial dimensions\n",
-    "# Also assume that 'da' has the original latitude and longitude data\n",
-    "latitudes = da.coords[\"latitude\"].values\n",
-    "longitudes = da.coords[\"longitude\"].values\n",
+    "da = dataset[\"Rrs\"]\n",
     "\n",
-    "# Filtering condition based on wavelength values\n",
     "filter_condition = (\n",
     "    (da.sel(wavelength=650) > da.sel(wavelength=620))\n",
     "    & (da.sel(wavelength=701) > da.sel(wavelength=681))\n",
     "    & (da.sel(wavelength=701) > da.sel(wavelength=450))\n",
     ")\n",
+    "extent = [-95, -85, 27, 33]\n",
+    "colors = [\"#e41a1c\", \"#377eb8\", \"#4daf4a\", \"#f781bf\", \"#a65628\", \"#984ea3\"]\n",
     "\n",
-    "# Apply the filtering condition to the SAM classification results\n",
-    "filtered_best_match_full = np.where(filter_condition, best_match_full, np.nan)\n",
-    "\n",
-    "# Create a custom discrete color map\n",
-    "cmap = mcolors.ListedColormap(\n",
-    "    [\"#e41a1c\", \"#377eb8\", \"#4daf4a\", \"#f781bf\", \"#a65628\", \"#984ea3\"]\n",
-    ")\n",
-    "bounds = np.arange(-0.5, n_clusters, 1)\n",
-    "norm = mcolors.BoundaryNorm(bounds, cmap.N)\n",
-    "\n",
-    "# Create a figure and axis with the correct map projection\n",
-    "fig, ax = plt.subplots(figsize=(10, 8), subplot_kw={\"projection\": ccrs.PlateCarree()})\n",
-    "\n",
-    "# Plot the filtered SAM classification results\n",
-    "im = ax.pcolormesh(\n",
-    "    longitudes,\n",
-    "    latitudes,\n",
-    "    filtered_best_match_full,\n",
-    "    cmap=cmap,\n",
-    "    norm=norm,\n",
-    "    transform=ccrs.PlateCarree(),\n",
-    ")\n",
-    "\n",
-    "# Add geographic features for context\n",
-    "ax.add_feature(cfeature.COASTLINE)\n",
-    "ax.add_feature(cfeature.BORDERS, linestyle=\":\")\n",
-    "ax.add_feature(cfeature.STATES, linestyle=\"--\")\n",
-    "\n",
-    "# Add gridlines\n",
-    "ax.gridlines(draw_labels=True)\n",
-    "\n",
-    "# Set the extent to zoom in to the specified region\n",
-    "ax.set_extent([-95, -85, 27, 33], crs=ccrs.PlateCarree())\n",
-    "\n",
-    "# Add color bar with labels\n",
-    "cbar = plt.colorbar(\n",
-    "    im,\n",
-    "    ax=ax,\n",
-    "    orientation=\"vertical\",\n",
-    "    pad=0.02,\n",
-    "    fraction=0.05,\n",
-    "    ticks=np.arange(n_clusters),\n",
-    ")\n",
-    "cbar.ax.set_yticklabels([f\"Class {i+1}\" for i in range(n_clusters)])\n",
-    "cbar.set_label(\"Water Types\", rotation=270, labelpad=20)\n",
-    "\n",
-    "# Add title\n",
-    "ax.set_title(\"PCA + Spectral Angle Mapper (SAM) Water Type Classification\", fontsize=16)\n",
-    "\n",
-    "# Show the plot\n",
-    "plt.show()"
+    "cluster_labels, latitudes, longitudes = apply_kmeans(\n",
+    "    da, n_clusters=6, filter_condition=filter_condition, extent=extent, colors=colors\n",
+    ")"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Spectral Angle Mapper (SAM)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import numpy as np\n",
-    "import pandas as pd\n",
-    "import matplotlib.pyplot as plt\n",
-    "import matplotlib.colors as mcolors\n",
-    "from scipy.interpolate import interp1d"
+    "## Principal Component Analysis (PCA)"
    ]
   },
   {
@@ -916,10 +251,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Load the netCDF data\n",
-    "file_path = \"data/PACE_OCI.20240730T181157.L2.OC_AOP.V2_0.NRT.nc\"\n",
-    "dataset = hypercoast.read_pace(file_path)\n",
-    "da = dataset[\"Rrs\"]  # Assuming 'Rrs' contains the reflectance data"
+    "pca_data = apply_pca(dataset, n_components=3, x_component=0, y_component=1)"
    ]
   },
   {
@@ -928,73 +260,16 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Extract PACE wavelengths\n",
-    "pace_wavelengths = da[\"wavelength\"].values\n",
-    "print(pace_wavelengths)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# Function to load and resample a single CSV spectral library file\n",
-    "def load_and_resample_spectral_library(csv_path, target_wavelengths):\n",
-    "    df = pd.read_csv(csv_path)\n",
-    "    original_wavelengths = df.iloc[:, 0].values  # First column is wavelength\n",
-    "    spectra_values = df.iloc[:, 1].values  # Second column is spectral values\n",
-    "\n",
-    "    # Interpolation function\n",
-    "    interp_func = interp1d(\n",
-    "        original_wavelengths, spectra_values, kind=\"linear\", fill_value=\"extrapolate\"\n",
-    "    )\n",
-    "\n",
-    "    # Resample to the target (PACE) wavelengths\n",
-    "    resampled_spectra = interp_func(target_wavelengths)\n",
-    "\n",
-    "    return resampled_spectra"
+    "pca_data = apply_pca(dataset, n_components=3, x_component=1, y_component=2)"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Download the SAM spectral library:"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "url = \"https://github.com/opengeos/datasets/releases/download/hypercoast/SAM_spectral_library.zip\"\n",
-    "hypercoast.download_file(url)"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# Load and resample all 6 endmembers\n",
-    "endmember_paths = [\n",
-    "    \"./SAM_spectral_library/Dataset_1.csv\",\n",
-    "    \"./SAM_spectral_library/Dataset_2.csv\",\n",
-    "    \"./SAM_spectral_library/Dataset_3.csv\",\n",
-    "    \"./SAM_spectral_library/Dataset_4.csv\",\n",
-    "    \"./SAM_spectral_library/Dataset_5.csv\",\n",
-    "    \"./SAM_spectral_library/Dataset_6.csv\",\n",
-    "]\n",
+    "## Spectral Angle Mapper (SAM)\n",
     "\n",
-    "endmembers = np.array(\n",
-    "    [\n",
-    "        load_and_resample_spectral_library(path, pace_wavelengths)\n",
-    "        for path in endmember_paths\n",
-    "    ]\n",
-    ")"
+    "### Apply SAM to the whole image"
    ]
   },
   {
@@ -1003,44 +278,18 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "print(endmembers)"
+    "data, latitudes, longitudes = apply_sam(\n",
+    "    dataset,\n",
+    "    n_components=3,\n",
+    "    n_clusters=6,\n",
+    ")"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "# Plot sample spectra from the CSV files and their resampled versions\n",
-    "def plot_sample_spectra(csv_paths, pace_wavelengths):\n",
-    "    plt.figure(figsize=(14, 8))\n",
-    "\n",
-    "    for i, csv_path in enumerate(csv_paths):\n",
-    "        df = pd.read_csv(csv_path)\n",
-    "        original_wavelengths = df.iloc[:, 0].values\n",
-    "        spectra_values = df.iloc[:, 1].values\n",
-    "        resampled_spectra = load_and_resample_spectral_library(\n",
-    "            csv_path, pace_wavelengths\n",
-    "        )\n",
-    "\n",
-    "        plt.plot(\n",
-    "            original_wavelengths,\n",
-    "            spectra_values,\n",
-    "            label=f\"Original Spectra {i+1}\",\n",
-    "            linestyle=\"--\",\n",
-    "        )\n",
-    "        plt.plot(pace_wavelengths, resampled_spectra, label=f\"Resampled Spectra {i+1}\")\n",
-    "\n",
-    "    plt.xlabel(\"Wavelength (nm)\")\n",
-    "    plt.ylabel(\"Spectral Reflectance\")\n",
-    "    plt.title(\"Comparison of Original and Resampled Spectra\")\n",
-    "    plt.legend()\n",
-    "    plt.grid(True)\n",
-    "    plt.show()\n",
-    "\n",
-    "\n",
-    "plot_sample_spectra(endmember_paths, pace_wavelengths)"
+    "### Apply SAM to selected pixels"
    ]
   },
   {
@@ -1049,23 +298,22 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Function to calculate spectral angle\n",
-    "def spectral_angle_mapper(pixel, reference):\n",
-    "    norm_pixel = np.linalg.norm(pixel)\n",
-    "    norm_reference = np.linalg.norm(reference)\n",
-    "    cos_theta = np.dot(pixel, reference) / (norm_pixel * norm_reference)\n",
-    "    angle = np.arccos(np.clip(cos_theta, -1, 1))\n",
-    "    return angle"
+    "extent = [-95, -85, 27, 33]\n",
+    "colors = [\"#377eb8\", \"#ff7f00\", \"#4daf4a\", \"#f781bf\", \"#a65628\", \"#984ea3\"]\n",
+    "data, latitudes, longitudes = apply_sam(\n",
+    "    dataset,\n",
+    "    n_components=3,\n",
+    "    n_clusters=6,\n",
+    "    extent=extent,\n",
+    "    colors=colors,\n",
+    ")"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "# Reshape data to (n_pixels, n_bands)\n",
-    "reshaped_data = da.values.reshape(-1, da.shape[-1])"
+    "### Apply SAM with a filtering condition"
    ]
   },
   {
@@ -1074,22 +322,31 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Apply SAM for each pixel and each endmember\n",
-    "angles = np.zeros((reshaped_data.shape[0], endmembers.shape[0]))\n",
+    "da = dataset[\"Rrs\"]\n",
     "\n",
-    "for i in range(reshaped_data.shape[0]):\n",
-    "    for j in range(endmembers.shape[0]):\n",
-    "        angles[i, j] = spectral_angle_mapper(reshaped_data[i, :], endmembers[j, :])"
+    "filter_condition = (\n",
+    "    (da.sel(wavelength=650) > da.sel(wavelength=620))\n",
+    "    & (da.sel(wavelength=701) > da.sel(wavelength=681))\n",
+    "    & (da.sel(wavelength=701) > da.sel(wavelength=450))\n",
+    ")\n",
+    "extent = [-95, -85, 27, 33]\n",
+    "colors = [\"#e41a1c\", \"#377eb8\", \"#4daf4a\", \"#f781bf\", \"#a65628\", \"#984ea3\"]\n",
+    "\n",
+    "data, latitudes, longitudes = apply_sam(\n",
+    "    dataset,\n",
+    "    n_components=3,\n",
+    "    n_clusters=6,\n",
+    "    filter_condition=filter_condition,\n",
+    "    extent=extent,\n",
+    "    colors=colors,\n",
+    ")"
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "metadata": {},
-   "outputs": [],
    "source": [
-    "# Find the minimum angle (best match) for each pixel\n",
-    "best_match = np.argmin(angles, axis=1)"
+    "### Use spectral library"
    ]
   },
   {
@@ -1098,8 +355,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Reshape best_match back to the original spatial dimensions\n",
-    "best_match = best_match.reshape(da.shape[:-1])"
+    "filepath = \"data/PACE_OCI.20240730T181157.L2.OC_AOP.V2_0.NRT.nc\"\n",
+    "dataset = hypercoast.read_pace(filepath)\n",
+    "url = \"https://github.com/opengeos/datasets/releases/download/hypercoast/SAM_spectral_library.zip\"\n",
+    "hypercoast.download_file(url)\n",
+    "spectral_library = \"./SAM_spectral_library/*.csv\""
    ]
   },
   {
@@ -1108,63 +368,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Assume 'best_match' contains the SAM classification results reshaped to the original spatial dimensions\n",
-    "# Also assume that 'da' has the original latitude and longitude data\n",
-    "latitudes = da.coords[\"latitude\"].values\n",
-    "longitudes = da.coords[\"longitude\"].values\n",
-    "\n",
-    "# Create a custom discrete color map\n",
-    "cmap = mcolors.ListedColormap(\n",
-    "    [\"#377eb8\", \"#e41a1c\", \"#4daf4a\", \"#f781bf\", \"#a65628\", \"#984ea3\"]\n",
-    ")\n",
-    "bounds = np.arange(-0.5, len(endmembers), 1)\n",
-    "norm = mcolors.BoundaryNorm(bounds, cmap.N)\n",
-    "\n",
-    "# Create a figure and axis with the correct map projection\n",
-    "fig, ax = plt.subplots(figsize=(12, 10), subplot_kw={\"projection\": ccrs.PlateCarree()})\n",
-    "\n",
-    "# Plot the SAM classification results\n",
-    "im = ax.pcolormesh(\n",
-    "    longitudes,\n",
-    "    latitudes,\n",
-    "    best_match,\n",
-    "    cmap=cmap,\n",
-    "    norm=norm,\n",
-    "    transform=ccrs.PlateCarree(),\n",
-    ")\n",
-    "\n",
-    "# Add geographic features for context\n",
-    "ax.add_feature(cfeature.COASTLINE)\n",
-    "ax.add_feature(cfeature.BORDERS, linestyle=\":\")\n",
-    "ax.add_feature(cfeature.STATES, linestyle=\"--\")\n",
-    "\n",
-    "# Adding axis labels\n",
-    "ax.set_xlabel(\"Longitude\")\n",
-    "ax.set_ylabel(\"Latitude\")\n",
-    "\n",
-    "# Adding a title\n",
-    "ax.set_title(\"Spectral Angle Mapper (SAM) Water Type Classification\", fontsize=16)\n",
-    "\n",
-    "# Adding a color bar with discrete values\n",
-    "cbar = plt.colorbar(\n",
-    "    im,\n",
-    "    ax=ax,\n",
-    "    orientation=\"vertical\",\n",
-    "    pad=0.02,\n",
-    "    fraction=0.05,\n",
-    "    ticks=np.arange(len(endmembers)),\n",
-    ")\n",
-    "cbar.ax.set_yticklabels([f\"Class {i+1}\" for i in range(len(endmembers))])\n",
-    "cbar.set_label(\"Water Types\", rotation=270, labelpad=20)\n",
-    "\n",
-    "# Adding gridlines\n",
-    "ax.gridlines(draw_labels=True, linestyle=\"--\", linewidth=0.5)\n",
-    "\n",
-    "# Set the extent to zoom in to the specified region (adjust as needed)\n",
-    "ax.set_extent([-95, -85, 27, 33], crs=ccrs.PlateCarree())\n",
-    "\n",
-    "# Show the plot\n",
-    "plt.show()"
+    "extent = [-95, -85, 27, 33]\n",
+    "data, latitudes, longitudes = apply_sam_spectral(\n",
+    "    dataset,\n",
+    "    spectral_library=spectral_library,\n",
+    "    extent=extent,\n",
+    ")"
    ]
   },
   {
@@ -1173,80 +382,19 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Assume 'best_match' contains the SAM classification results reshaped to the original spatial dimensions\n",
-    "# Also assume that 'da' has the original latitude and longitude data\n",
-    "latitudes = da.coords[\"latitude\"].values\n",
-    "longitudes = da.coords[\"longitude\"].values\n",
-    "\n",
-    "# Extract specific wavelengths for the conditions\n",
-    "wavelength_450 = da.sel(wavelength=450).values\n",
-    "wavelength_620 = da.sel(wavelength=620).values\n",
-    "wavelength_650 = da.sel(wavelength=650).values\n",
-    "wavelength_681 = da.sel(wavelength=681).values\n",
-    "wavelength_701 = da.sel(wavelength=701).values\n",
-    "\n",
-    "# Apply the condition to filter the pixels\n",
-    "condition = (\n",
-    "    (wavelength_650 > wavelength_620)\n",
-    "    & (wavelength_701 > wavelength_681)\n",
-    "    & (wavelength_701 > wavelength_450)\n",
-    ")\n",
-    "\n",
-    "# Filter the best_match data based on the condition\n",
-    "filtered_best_match = np.where(condition, best_match, np.nan)\n",
-    "\n",
-    "# Create a custom discrete color map\n",
-    "cmap = mcolors.ListedColormap(\n",
-    "    [\"#377eb8\", \"#e41a1c\", \"#4daf4a\", \"#f781bf\", \"#a65628\", \"#984ea3\"]\n",
-    ")\n",
-    "bounds = np.arange(-0.5, len(endmembers), 1)\n",
-    "norm = mcolors.BoundaryNorm(bounds, cmap.N)\n",
-    "\n",
-    "# Create a figure and axis with the correct map projection\n",
-    "fig, ax = plt.subplots(figsize=(12, 10), subplot_kw={\"projection\": ccrs.PlateCarree()})\n",
-    "\n",
-    "# Plot the SAM classification results for the filtered pixels\n",
-    "im = ax.pcolormesh(\n",
-    "    longitudes,\n",
-    "    latitudes,\n",
-    "    filtered_best_match,\n",
-    "    cmap=cmap,\n",
-    "    norm=norm,\n",
-    "    transform=ccrs.PlateCarree(),\n",
-    ")\n",
-    "\n",
-    "# Add geographic features for context\n",
-    "ax.add_feature(cfeature.COASTLINE)\n",
-    "ax.add_feature(cfeature.BORDERS, linestyle=\":\")\n",
-    "ax.add_feature(cfeature.STATES, linestyle=\"--\")\n",
-    "\n",
-    "# Adding axis labels\n",
-    "ax.set_xlabel(\"Longitude\")\n",
-    "ax.set_ylabel(\"Latitude\")\n",
-    "\n",
-    "# Adding a title\n",
-    "ax.set_title(\"SAM Water Type Classification\", fontsize=16)\n",
-    "\n",
-    "# Adding a color bar with discrete values\n",
-    "cbar = plt.colorbar(\n",
-    "    im,\n",
-    "    ax=ax,\n",
-    "    orientation=\"vertical\",\n",
-    "    pad=0.02,\n",
-    "    fraction=0.05,\n",
-    "    ticks=np.arange(len(endmembers)),\n",
+    "da = dataset[\"Rrs\"]\n",
+    "extent = [-95, -85, 27, 33]\n",
+    "filter_condition = (\n",
+    "    (da.sel(wavelength=650) > da.sel(wavelength=620))\n",
+    "    & (da.sel(wavelength=701) > da.sel(wavelength=681))\n",
+    "    & (da.sel(wavelength=701) > da.sel(wavelength=450))\n",
     ")\n",
-    "cbar.ax.set_yticklabels([f\"Class {i+1}\" for i in range(len(endmembers))])\n",
-    "cbar.set_label(\"Water Types\", rotation=270, labelpad=20)\n",
-    "\n",
-    "# Adding gridlines\n",
-    "ax.gridlines(draw_labels=True, linestyle=\"--\", linewidth=0.5)\n",
-    "\n",
-    "# Set the extent to zoom in to the specified region (adjust as needed)\n",
-    "ax.set_extent([-95, -85, 27, 33], crs=ccrs.PlateCarree())\n",
-    "\n",
-    "# Show the plot\n",
-    "plt.show()"
+    "data, latitudes, longitudes = apply_sam_spectral(\n",
+    "    da,\n",
+    "    spectral_library=spectral_library,\n",
+    "    filter_condition=filter_condition,\n",
+    "    extent=extent,\n",
+    ")"
    ]
   }
  ],
@@ -1255,7 +403,7 @@
    "provenance": []
   },
   "kernelspec": {
-   "display_name": "Hyper",
+   "display_name": "hyper",
    "language": "python",
    "name": "python3"
   },