{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CryoTEMPO: Assessing Ice and Sea Level Dynamics near Petermann Glacier, Greenland"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction\n",
    "\n",
    "In this example notebook, we use the [CryoTEMPO](http://cryosat.mssl.ucl.ac.uk/tempo/) Land Ice, Sea Ice and Polar Oceans data stored in Specklia to investigate the dynamics of Petermann Glacier and its surrounding environment over 12 years.\n",
    "\n",
    "Petermann Glacier, one of Greenland's largest marine-terminating glaciers, has experienced significant thinning in recent decades, leading to substantial loss of its floating ice shelf through calving events. These changes hold implications for sea level rise.\n",
    "\n",
    "Petermann Glacier discharges into the Nares Strait, a critical pathway for sea ice exiting the Arctic. The glacier's ice shelf is seasonally shielded by sea ice in this pathway. However, with the thinning of sea ice, it becomes more susceptible to faster floes, potentially increasing glacier calving rates.\n",
    "\n",
    "This notebook aims to compare the long-term trends of sea-level fluctuations, sea ice thickness, and land ice thickness in and around Petermann Glacier using CryoSat's [CryoTEMPO](http://cryosat.mssl.ucl.ac.uk/tempo/) products.\n",
    "\n",
    "You can view the code below or [click here](https://colab.research.google.com/github/earthwave/specklia_demo_notebooks/blob/main/cryotempo_petermann_glacier.ipynb/) to run it yourself in Google Colab!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Environment Setup\n",
    "To run this notebook, you will need to make sure that the folllowing packages are installed in your python environment (all can be installed via pip/conda):\n",
    "- matplotlib\n",
    "- geopandas\n",
    "- contextily\n",
    "- ipykernel\n",
    "- shapely\n",
    "- specklia\n",
    "\n",
    "If you are using the Google Colab environment, these packages will be installed in the next cell. Please note this step may take a few minutes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "if 'google.colab' in sys.modules:\n",
    "    %pip install rasterio --no-binary rasterio\n",
    "    %pip install specklia\n",
    "    %pip install matplotlib\n",
    "    %pip install geopandas\n",
    "    %pip install contextily\n",
    "    %pip install shapely\n",
    "    %pip install python-dotenv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbsphinx": "hidden"
   },
   "outputs": [],
   "source": [
    "# fix an issue that can sometimes occur with rasterio using the wrong version of proj\n",
    "import os\n",
    "import pyproj\n",
    "os.environ['PROJ_LIB'] = pyproj.datadir.get_data_dir()\n",
    "\n",
    "from datetime import datetime\n",
    "import math\n",
    "import pprint\n",
    "from time import perf_counter\n",
    "from typing import Tuple\n",
    "\n",
    "import contextily as ctx\n",
    "from dotenv import load_dotenv\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from numpy.typing import NDArray\n",
    "import shapely\n",
    "from specklia import Specklia\n",
    "\n",
    "# load a demonstration API key from a .env file.\n",
    "load_dotenv()\n",
    "\n",
    "# Specify polygon colors to match another Earthwave product, https://cs2eo.org.\n",
    "data_coverage_color = np.array([52, 211, 153])/255\n",
    "geo_filter_colour = np.array([96, 165, 250])/255\n",
    "\n",
    "# local EPSG projection code for plotting, we're considering data in the northern hemisphere so will use the\n",
    "# Northern Polar Stereographic projection\n",
    "nps_epsg_code = 3413\n",
    "\n",
    "# To run this code yourself, first generate your own key using https://specklia.earthwave.co.uk.\n",
    "if 'DEMO_API_KEY' in os.environ:\n",
    "    client = Specklia(os.environ['DEMO_API_KEY'])\n",
    "else:\n",
    "    user_api_key = input('Please generate your own key using https://specklia.earthwave.co.uk/ApiKeys and paste it here:')\n",
    "    client = Specklia(user_api_key)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot dataset spatial coverage\n",
    "Let's plot the geospatial coverage of the CryoTEMPO Polar Oceans, Sea Ice and Land Ice datasets, using [Contextily](https://contextily.readthedocs.io/en/latest/) to add a satellite background map.\n",
    "All background satellite imagery tiles are provided by the [Esri World Imagery Service](https://www.esriuk.com/en-gb/content/products?esri-world-imagery-service) via Contextily."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "available_datasets = client.list_datasets()\n",
    "\n",
    "cryotempo_datasets = {}\n",
    "dataset_names = ['CryoTEMPO Polar Ocean Thematic Product Baseline B',\n",
    "                 'CryoTEMPO Sea Ice Thematic Product Baseline B',\n",
    "                 'CryoTEMPO Land Ice Thematic Product Baseline B']\n",
    "\n",
    "for ds_name in dataset_names:\n",
    "    cryotempo_datasets[ds_name] = available_datasets[\n",
    "        available_datasets['dataset_name'] == ds_name].iloc[0]\n",
    "\n",
    "    print(f\"{ds_name} contains data timestamped between \\n\"\n",
    "        f\"{cryotempo_datasets[ds_name]['min_timestamp']} \"\n",
    "        f\"and {cryotempo_datasets[ds_name]['max_timestamp']}\\n\")\n",
    "\n",
    "    print(f\"{ds_name} has the following columns:\")\n",
    "    pprint.PrettyPrinter(indent=2, width=120).pprint(cryotempo_datasets[ds_name]['columns'])\n",
    "    print(\"\\n\\n\")\n",
    "\n",
    "    desired_dataset_spatial_coverage = gpd.GeoDataFrame(\n",
    "        geometry=[cryotempo_datasets[ds_name]['epsg4326_coverage']], crs=4326)\n",
    "    \n",
    "    # we create two plots, one for each hemisphere, to minimise distortion\n",
    "    # and illustrate that multiple regions may be present within the same dataset\n",
    "    for crs, hemisphere_name, hemisphere_bounding_box in [\n",
    "            (3413, 'northern', (-180, 0, 180, 90)),\n",
    "            (3031, 'southern', (-180, -90, 180, 0))\n",
    "            ]:\n",
    "\n",
    "        cropped_data = desired_dataset_spatial_coverage.clip(hemisphere_bounding_box)\n",
    "\n",
    "        if len(cropped_data) > 0:\n",
    "            transformed_cropped_data = cropped_data.to_crs(crs)\n",
    "            ax = transformed_cropped_data.plot(\n",
    "                figsize=(10, 10), alpha=0.5, color=data_coverage_color, edgecolor=data_coverage_color)\n",
    "\n",
    "            ax.set_xlabel('x (m)')\n",
    "            ax.set_ylabel('y (m)')\n",
    "            ax.set_title(f\"{ds_name} spatial coverage, {hemisphere_name} hemisphere (EPSG {crs})\")\n",
    "\n",
    "            ctx.add_basemap(ax, source=ctx.providers.Esri.WorldImagery, crs=crs, attribution=False, zoom=2)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define area of interest\n",
    "\n",
    "Lets define a rough boundary encompassing the terminus of Petermann Glacier and the adjacent ocean.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "petermann_extent = shapely.Polygon((\n",
    "    (-64.53, 80.99), (-63.87, 81.16), (-62.84, 81.19), (-61.9, 81.02), (-61.3, 80.65),\n",
    "    (-59., 80.2), (-58., 80.4), (-60., 80.7), (-61.73, 81.29), (-61.3, 81.51),\n",
    "    (-61.71, 81.78), (-62.86, 81.97), (-67.4, 81.42), (-64.53, 80.99)))\n",
    "petermann_extent_gdf = gpd.GeoSeries(petermann_extent, crs=4326)\n",
    "\n",
    "ax = petermann_extent_gdf.to_crs(nps_epsg_code).plot(\n",
    "    figsize=(10, 10), alpha=0.5, color=geo_filter_colour, edgecolor=geo_filter_colour)\n",
    "\n",
    "ax.set_xlabel('x (meters)')\n",
    "ax.set_ylabel('y (meters)')\n",
    "ax.set_title(\"Query extent covering Peterman Glacier and Surrounding Ocean\")\n",
    "ctx.add_basemap(ax, source=ctx.providers.Esri.WorldImagery, crs=nps_epsg_code, attribution=False, zoom=8)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Query data\n",
    "Now that we've defined our query extent, lets retrieve data from the CryoTEMPO Polar Ocean Thematic Product, CryoTEMPO Sea Ice Thematic Product, and CryoTEMPO Land Ice Thematic Product from Specklia.\n",
    "\n",
    "Whilst we do so, we'll apply some uncertainty filtering and filter the CryoTEMPO Land Ice data on surface type to retrieve measurements from grounded (surface_type=1) and floating ice (surface_type=2) only. Note that the surface type definitions can be found in the [CryoTEMPO Product Handbook](http://www.cpom.ucl.ac.uk/cryotempo/pdf_viewer.php?theme=landice)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cryotempo_data = {}\n",
    "sources = {}\n",
    "query_parameters = {\n",
    "    'CryoTEMPO Polar Ocean Thematic Product Baseline B': {\n",
    "        'filters': [{'column': 'sea_level_anomaly_uncertainty', 'operator': '<=', 'threshold': .05}],\n",
    "        'columns': ['sea_level_anomaly_filtered', 'timestamp']},\n",
    "    'CryoTEMPO Sea Ice Thematic Product Baseline B': {\n",
    "        'filters': [{'column': 'sea_ice_freeboard_uncertainty', 'operator': '<=', 'threshold': .15}],\n",
    "        'columns': ['sea_ice_freeboard_filtered', 'timestamp']},\n",
    "    'CryoTEMPO Land Ice Thematic Product Baseline B': {\n",
    "        'filters': [{'column': 'uncertainty', 'operator': '<=', 'threshold': 1.5},\n",
    "                    {'column': 'surface_type', 'operator': '>=', 'threshold': 1.},\n",
    "                    {'column': 'surface_type', 'operator': '<=', 'threshold': 2.}],\n",
    "        'columns': ['elevation', 'reference_dem', 'timestamp', 'surface_type']}}\n",
    "\n",
    "for ds_name in dataset_names:\n",
    "    dataset = available_datasets[available_datasets['dataset_name'] == ds_name].iloc[0]\n",
    "    query_start_time = perf_counter()\n",
    "    cryotempo_data[ds_name], sources[ds_name] = client.query_dataset(\n",
    "        dataset_id=dataset['dataset_id'],\n",
    "        epsg4326_polygon=petermann_extent,\n",
    "        min_datetime=datetime(2011, 1, 1),\n",
    "        max_datetime=datetime(2023, 12, 1),\n",
    "        columns_to_return=query_parameters[ds_name]['columns'],\n",
    "        additional_filters=query_parameters[ds_name]['filters'])\n",
    "\n",
    "    print(f'Query took {perf_counter()-query_start_time:.2f} seconds to complete.')\n",
    "    print(f'{ds_name} Query complete, {len(cryotempo_data[ds_name])} points returned, '\n",
    "            f'drawn from {len(sources[ds_name])} original sources.')\n",
    "    print(f'Columns within the data: {cryotempo_data[ds_name].columns}\\n\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analyse short and long-term trends\n",
    "\n",
    "We next define some functions to help us disseminate the point data we've retrieved and analyse short and long term trends in the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_linear_trend(point_data_gdf: gpd.GeoDataFrame, column_name: str) -> Tuple[NDArray, NDArray, str]:\n",
    "    \"\"\"\n",
    "    Find the line of best fit for the point data using least squares linear regression.\n",
    "\n",
    "    Format an associated label for the trend line in appropriate units for readability.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    point_data_gdf : gpd.GeoDataFrame\n",
    "        Point data for data set.\n",
    "    column_name : str\n",
    "        Column in point data geodataframe to apply least squares method to.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    Tuple[NDArray, NDArray, str]\n",
    "        The X-values, Y-values, and label for plotting the trend line.\n",
    "    \"\"\"\n",
    "    # remove nan values to avoid propagation in regression\n",
    "    filtered_point_data_gdf = point_data_gdf.loc[point_data_gdf[column_name].notna()].copy()\n",
    "\n",
    "    # trend line\n",
    "    m, c = np.linalg.lstsq(\n",
    "        np.vstack([filtered_point_data_gdf.timestamp.to_numpy(), np.ones(len(filtered_point_data_gdf))]).T,\n",
    "        filtered_point_data_gdf[column_name].to_numpy(),\n",
    "        rcond=None)[0]\n",
    "    # generate x coordinates (timestamps in seconds) in 30 day intervals\n",
    "    x_coords = np.arange(filtered_point_data_gdf.timestamp.min(), filtered_point_data_gdf.timestamp.max(), 60*60*24*30)\n",
    "    y_trend_line = m * x_coords + c\n",
    "    x_datetime = [datetime.fromtimestamp(x) for x in x_coords]\n",
    "\n",
    "    # format trend units to make them more readable\n",
    "    meter_per_year_trend = m * 60 * 60 * 24 * 365\n",
    "    order_of_magnitude = math.floor(math.log(abs(meter_per_year_trend), 10))\n",
    "    if order_of_magnitude <= -3:\n",
    "        unit = 'mm'\n",
    "        rounded_scaled_trend = round(meter_per_year_trend / 10**order_of_magnitude, 2)\n",
    "    elif order_of_magnitude == -2:\n",
    "        unit = 'cm'\n",
    "        rounded_scaled_trend = round(meter_per_year_trend / 10**order_of_magnitude, 2)\n",
    "    else:\n",
    "        unit = 'm'\n",
    "        rounded_scaled_trend = round(meter_per_year_trend, 2)\n",
    "\n",
    "    label = f'trend: {rounded_scaled_trend} {unit}/y'\n",
    "    return x_datetime, y_trend_line, label\n",
    "\n",
    "\n",
    "def calculate_monthly_average_time_series(point_data_gdf: gpd.GeoDataFrame, column_name: str) -> gpd.GeoDataFrame:\n",
    "    \"\"\"\n",
    "    Calculate monthly average of point data for full timespan.\n",
    "\n",
    "    Note that centre date of each month has been set to the 15th for simplicity.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    point_data_gdf : gpd.GeoDataFrame\n",
    "        Point data for data set.\n",
    "    column_name : str\n",
    "        Column in point data geodataframe to calculate monthly average of.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    gpd.GeoDataFrame\n",
    "        Contains monthly average of selected column and centre date of each month.\n",
    "    \"\"\"\n",
    "    point_data_gdf['centre_month'] = [\n",
    "        datetime(datetime.fromtimestamp(t).year, datetime.fromtimestamp(t).month, 15) for t in point_data_gdf.timestamp]\n",
    "    filtered_point_data_gdf = point_data_gdf.loc[point_data_gdf[column_name].notna()].copy()\n",
    "    monthly_average_df = filtered_point_data_gdf.groupby(by='centre_month').agg({column_name: 'mean'}).reset_index()\n",
    "    return monthly_average_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using these functions we've defined we can calculate and plot changes in sea level, sea ice thickness, ice shelf elevation and land ice elevation over time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(figsize=(15, 5 * (len(cryotempo_data) + 1)),\n",
    "                         nrows=len(cryotempo_data) + 1)\n",
    "\n",
    "for i, (ds_name, gdf) in enumerate(cryotempo_data.items()):\n",
    "    # define our columns of interest in the different datasets\n",
    "    if ds_name == 'CryoTEMPO Polar Ocean Thematic Product Baseline B':\n",
    "        column_name = 'sea_level_anomaly_filtered'\n",
    "    elif ds_name == 'CryoTEMPO Sea Ice Thematic Product Baseline B':\n",
    "        column_name = 'sea_ice_freeboard_filtered'\n",
    "    elif ds_name == 'CryoTEMPO Land Ice Thematic Product Baseline B':\n",
    "        column_name = 'elevation_difference'\n",
    "        gdf['elevation_difference'] = gdf['elevation'] - gdf['reference_dem']\n",
    "    \n",
    "    # remove underscores and capitalise for label\n",
    "    y_label = column_name.replace(\"_\", \" \").title()\n",
    "\n",
    "    if ds_name in ['CryoTEMPO Polar Ocean Thematic Product Baseline B', 'CryoTEMPO Sea Ice Thematic Product Baseline B']:\n",
    "        # plot point data in background\n",
    "        axes[i].scatter([datetime.fromtimestamp(t) for t in gdf.timestamp],\n",
    "                        gdf[column_name], s=0.5, c='lightgrey', label='point data')\n",
    "\n",
    "        # calculate and plot monthly average\n",
    "        monthly_average_df = calculate_monthly_average_time_series(gdf, column_name)\n",
    "        axes[i].plot(monthly_average_df['centre_month'], monthly_average_df[column_name], label='monthly average')\n",
    "\n",
    "        # calcualte and plot linear trend\n",
    "        x_datetime, y_trend_line, trend_label = calculate_linear_trend(gdf, column_name)\n",
    "        axes[i].plot(x_datetime, y_trend_line, label=trend_label, c='red', ls='--')\n",
    "\n",
    "        axes[i].set_title(f'Time series of {y_label} from the {ds_name}')\n",
    "        ax.set_ylabel(y_label + ' [metres]')\n",
    "\n",
    "    elif ds_name == 'CryoTEMPO Land Ice Thematic Product Baseline B':\n",
    "        # we need to separate the floating ice from the grounded ice for the time series calculation\n",
    "        for j, (surface_type_name, surface_type_id) in enumerate({'grounded ice': 1, 'floating ice': 2}.items()):\n",
    "            surface_type_gdf = gdf.loc[gdf['surface_type'] == surface_type_id].copy()\n",
    "            # plot point data in background\n",
    "            axes[i + j].scatter([datetime.fromtimestamp(t) for t in surface_type_gdf.timestamp],\n",
    "                            surface_type_gdf[column_name], s=0.5, c='lightgrey', label='point data')\n",
    "\n",
    "            # calculate and plot monthly average\n",
    "            monthly_average_df = calculate_monthly_average_time_series(surface_type_gdf, column_name)\n",
    "            axes[i + j].plot(monthly_average_df['centre_month'], monthly_average_df[column_name],\n",
    "                             label='monthly average')\n",
    "\n",
    "            # calcualte and plot linear trend\n",
    "            x_datetime, y_trend_line, trend_label = calculate_linear_trend(surface_type_gdf, column_name)\n",
    "            axes[i + j].plot(x_datetime, y_trend_line, label=trend_label, c='red', ls='--')\n",
    "            axes[i + j].set_ylim(-15, 15)\n",
    "            axes[i + j].set_title(f'Time series of {y_label} over {surface_type_name} from the {ds_name}')\n",
    "            ax.set_ylabel(y_label + ' [metres]')\n",
    "\n",
    "# format axes\n",
    "for ax in axes:\n",
    "    ax.legend(loc='upper right')\n",
    "    ax.axhline(y=0, ls='--', c='black')\n",
    "    ax.set_xlabel('Date')\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The monthly time series above allows us to see short term variations in each variable. Most notably, we see clear seasonality in the thickening of grounded ice (bottom left) and floating ice (bottom right) in winter months and thinning in summer months. Sea ice freeboard is only measured between October and April each year, however we see how the sea ice generally thickens and then thins through that time.\n",
    "\n",
    "The long-term trend lines indicate a local sea level rise of 7.6 milimeter per year, with sea ice thinning by -1.6 milimeter per year, the Petermann Glacier ice shelf thinning by -0.3 metres per year, and grounded ice by 0.9 metres per year."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Map trends spatially\n",
    "\n",
    "It can also be useful to visualise the rate of change of each of these variables spatially.\n",
    "\n",
    "Below we define a function to calculate the linear trend through time at 2 kilometre resolution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_rate_of_change_map(\n",
    "        point_data_gdf: gpd.GeoDataFrame, column_name: str, resolution: int, epsg_code: int) -> gpd.GeoDataFrame:\n",
    "    \"\"\"\n",
    "    Calculate the rate of change per pixel in years.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    point_data_gdf : gpd.GeoDataFrame\n",
    "        Point data for data set.\n",
    "    column_name : str\n",
    "        Column name of point data to calculate rate of change over.\n",
    "    resolution : int\n",
    "        Target resolution of rate of change map.\n",
    "    epsg_code : int\n",
    "        Projection code with units of measurement in metres. \n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    gpd.GeoDataFrame\n",
    "        Rate of change in years per pixel of defined resolution.\n",
    "    \"\"\"\n",
    "    # reproject data so that x and y coordinates are in metres\n",
    "    df_reprojected = point_data_gdf.to_crs(epsg_code)\n",
    "\n",
    "    # group points per pixel\n",
    "    df_reprojected['x_centre'] = [math.floor(x / resolution) * resolution + resolution / 2\n",
    "                                  for x in df_reprojected.geometry.x]\n",
    "    df_reprojected['y_centre'] = [math.floor(y / resolution) * resolution + resolution / 2\n",
    "                                  for y in df_reprojected.geometry.y]\n",
    "    pixel_groups = df_reprojected.groupby(by=['x_centre', 'y_centre'])\n",
    "\n",
    "    # apply least squares linear regression to data in each pixel\n",
    "    rate_of_change = []\n",
    "    pixel_geometries = []\n",
    "    for pixel_coordinate, pixel_group in pixel_groups:\n",
    "        pixel_group_df = pixel_group.sort_values(by='timestamp')\n",
    "        pixel_group_df = pixel_group_df.loc[pixel_group_df[column_name].notna()]\n",
    "        if pixel_group_df.empty:\n",
    "            continue\n",
    "        m, _ = np.linalg.lstsq(\n",
    "            np.vstack([pixel_group_df['timestamp'].to_numpy(), np.ones(len(pixel_group_df))]).T,\n",
    "            pixel_group_df[column_name].to_numpy(),\n",
    "            rcond=None)[0]\n",
    "\n",
    "        # m gives the rate of change in meters per second (as timestamp is seconds), lets crudely convert it to\n",
    "        # metres per years for a more meaningful value\n",
    "        rate_of_change_meters_per_year = m * 60 * 60 * 24 * 365\n",
    "\n",
    "        # append results to list\n",
    "        rate_of_change.append(rate_of_change_meters_per_year)\n",
    "        pixel_geometries.append(shapely.Polygon((\n",
    "            (pixel_coordinate[0] + resolution / 2, pixel_coordinate[1] + resolution / 2),\n",
    "            (pixel_coordinate[0] + resolution / 2, pixel_coordinate[1] - resolution / 2),\n",
    "            (pixel_coordinate[0] - resolution / 2, pixel_coordinate[1] - resolution / 2),\n",
    "            (pixel_coordinate[0] - resolution / 2, pixel_coordinate[1] + resolution / 2),\n",
    "            (pixel_coordinate[0] + resolution / 2, pixel_coordinate[1] + resolution / 2))))\n",
    "\n",
    "    # create geodataframe containing gridded rate of change values\n",
    "    rate_of_change_gdf = gpd.GeoDataFrame(\n",
    "        {f'Rate of {column_name.replace(\"_\", \" \").title()} Change [m/y]': rate_of_change},\n",
    "        geometry=pixel_geometries, crs=epsg_code)\n",
    "    return rate_of_change_gdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can calculate and plot these rate of change maps for each dataset where, most notably, we see the significant thinning of Peterman Glacier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# target resolution of rate of change map in metres\n",
    "resolution = 2000\n",
    "\n",
    "fig, axes = plt.subplots(figsize=(5 * len(cryotempo_data), 7.5), ncols=len(cryotempo_data))\n",
    "min_x, min_y, max_x, max_y = petermann_extent_gdf.to_crs(nps_epsg_code).total_bounds\n",
    "for i, (ds_name, gdf) in enumerate(cryotempo_data.items()):\n",
    "    # define columns to fit regression to\n",
    "    if ds_name == 'CryoTEMPO Polar Ocean Thematic Product Baseline B':\n",
    "        column_name = 'sea_level_anomaly_filtered'\n",
    "    elif ds_name == 'CryoTEMPO Sea Ice Thematic Product Baseline B':\n",
    "        column_name = 'sea_ice_freeboard_filtered'\n",
    "    elif ds_name == 'CryoTEMPO Land Ice Thematic Product Baseline B':\n",
    "        column_name = 'elevation'\n",
    "    \n",
    "    # calculate rate of change map\n",
    "    rate_of_change_gdf = calculate_rate_of_change_map(gdf, column_name, resolution, nps_epsg_code)\n",
    "\n",
    "    # plot rate of change map\n",
    "    rate_of_change_col = rate_of_change_gdf.columns[0]\n",
    "    std = rate_of_change_gdf[rate_of_change_col].std()\n",
    "\n",
    "    # define colorbar bounds and spacing for tick labels\n",
    "    order_of_magnitude = math.floor(math.log(std, 10))\n",
    "    vmin = round(-2*std, abs(order_of_magnitude))\n",
    "    vmax = round(2*std, abs(order_of_magnitude))\n",
    "    spacing = round((vmax - vmin) / 4, abs(order_of_magnitude))\n",
    "\n",
    "    rate_of_change_gdf.plot(\n",
    "        rate_of_change_col, ax=axes[i], legend=True, vmin=vmin, vmax=vmax, cmap='RdYlBu',\n",
    "        legend_kwds={'label': rate_of_change_col, 'orientation': 'horizontal', 'shrink': .8,\n",
    "                     'pad': .1, 'ticks': np.arange(vmin, vmax + spacing, spacing)})\n",
    "    \n",
    "    # add background map\n",
    "    axes[i].set_xlim(min_x - resolution, max_x + resolution)\n",
    "    axes[i].set_ylim(min_y - resolution, max_y + resolution)\n",
    "    ctx.add_basemap(axes[i], source=ctx.providers.Esri.WorldImagery, crs=nps_epsg_code,\n",
    "                    attribution=False, zoom=8)\n",
    "\n",
    "    # format axis\n",
    "    axes[i].ticklabel_format(scilimits=(0, 2))\n",
    "    axes[i].set_xlabel('x (meters)')\n",
    "    axes[i].set_ylabel('y (meters)')\n",
    "    axes[i].set_title(ds_name)\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "specklia",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}