{
"cells": [
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"# CryoTEMPO-EOLIS Gridded Product & RGI: Comparing Elevation Changes of Surging and Non-Surging Glaciers in Svalbard"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"### Introduction\n",
"This example retrieves [CryoTEMPO-EOLIS Gridded Product](https://cryotempo-eolis.org/gridded-product/) data over 12 years using the [RGI v7.0](https://www.glims.org/RGI/) glacier boundary definitions, and analyses the elevation changes for different surge-type glaciers across Svalbard.\n",
"\n",
"You can view the code below or [click here](https://colab.research.google.com/github/earthwave/specklia_demo_notebooks/blob/main/eolis_svalbard_timeseries.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",
"from time import perf_counter\n",
"from typing import List, Tuple\n",
"\n",
"import contextily as ctx\n",
"from dotenv import load_dotenv\n",
"import geopandas as gpd\n",
"from matplotlib import cm\n",
"import matplotlib.gridspec as gridspec\n",
"from matplotlib.patches import Patch\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from numpy.typing import NDArray\n",
"from pandas import Timedelta\n",
"import shapely\n",
"from specklia import Specklia\n",
"\n",
"# load a demonstration API key from a .env file.\n",
"load_dotenv()\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": [
"### Query RGI polygons\n",
"\n",
"Firstly, let's use Specklia to retrieve the RGI polygons for glaciers in Svalbard."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"svalbard_polygon = shapely.Polygon(([3.9, 79.88], [28.6, 80.97], [28.6, 77.03], [14.57, 75.51], [3.9, 79.88]))\n",
"\n",
"available_datasets = client.list_datasets()\n",
"rgiv7_dataset = available_datasets[available_datasets['dataset_name'] == 'RGIv7 Glaciers'].iloc[0]\n",
"rgi_data, rgi_data_sources = client.query_dataset(\n",
" dataset_id=rgiv7_dataset['dataset_id'],\n",
" epsg4326_polygon=svalbard_polygon,\n",
" min_datetime=rgiv7_dataset['min_timestamp'] - Timedelta(seconds=1),\n",
" max_datetime=rgiv7_dataset['max_timestamp'] + Timedelta(seconds=1))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The RGI dataset provides a __surge_type__ attribute categorising each glacier based on evidence of surging. Surging is typically short-lived events where a glacier moves many times its normal rate.\n",
"\n",
"From the [RGI v7.0 user guide](https://zenodo.org/records/8362857) we get the __surge_type__ definitions:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"surge_type_definitions = {0: 'No evidence', 1: 'Possible', 2: 'Probable', 3: 'Observed', 9: 'Not assigned'}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can now plot the glacier boundaries on a map and colour code them by surge-type to see where they are."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots(figsize=(10, 10))\n",
"legend_elements = []\n",
"glacier_polygons = {}\n",
"for i, (surge_id, surge_type) in enumerate(surge_type_definitions.items()):\n",
" # we use a simplified version of the RGI polygons, which has a reduced number of vertices\n",
" polygons_for_surge_type = gpd.GeoSeries(\n",
" rgi_data[rgi_data.surge_type == surge_id]['simplified_polygons_t200'], crs=rgi_data.crs)\n",
"\n",
" if not polygons_for_surge_type.empty:\n",
" glacier_polygons[surge_type] = polygons_for_surge_type\n",
"\n",
" color = cm.CMRmap((i + 1) / (len(surge_type_definitions)))\n",
" polygons_for_surge_type.to_crs(nps_epsg_code).plot(ax=ax, facecolor=color, edgecolor=color, alpha=0.5)\n",
" legend_elements.append(Patch(facecolor=color, edgecolor=color,\n",
" label=surge_type + f' ({len(polygons_for_surge_type)})',\n",
" alpha=0.5))\n",
"ax.set_xlabel('x [m]')\n",
"ax.set_ylabel('y [m]')\n",
"ax.legend(handles=legend_elements, loc='upper right', title='Surge type (count)')\n",
"ctx.add_basemap(ax, source=ctx.providers.Esri.WorldImagery, crs=nps_epsg_code, zoom=6)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Query data for each surge-type category\n",
"\n",
"Now that we have the glacier boundary definitions for each surge-type we can query the CryoTEMPO-EOLIS gridded dataset."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"svalbard_eolis_data = {}\n",
"svalbard_data_sources = {}\n",
"\n",
"dataset_name = 'CryoTEMPO-EOLIS Gridded Product'\n",
"gridded_product_dataset = available_datasets[\n",
" available_datasets['dataset_name'] == dataset_name].iloc[0]\n",
"\n",
"for surge_id, surge_type in surge_type_definitions.items():\n",
" if surge_type not in glacier_polygons:\n",
" print(f'No data to query as there are no glaciers of surge type: {surge_type}, skipping..')\n",
" continue\n",
" query_polygon = glacier_polygons[surge_type].unary_union\n",
"\n",
" query_start_time = perf_counter()\n",
" svalbard_eolis_data[surge_type], svalbard_data_sources[surge_type] = client.query_dataset(\n",
" dataset_id=gridded_product_dataset['dataset_id'],\n",
" epsg4326_polygon=query_polygon,\n",
" min_datetime=datetime(2011, 1, 1),\n",
" max_datetime=datetime(2023, 11, 1),\n",
" columns_to_return=['timestamp', 'elevation', 'uncertainty'])\n",
"\n",
" print(f'Query of glaciers of surge type: {surge_type} complete in '\n",
" f'{perf_counter()-query_start_time:.2f} seconds, '\n",
" f'{len(svalbard_eolis_data[surge_type])} points returned, '\n",
" f'drawn from {len(svalbard_data_sources[surge_type])} original sources.')"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"### Compare elevation change time series per surge-type\n",
"\n",
"From the Gridded Product [documentation](https://cryotempo-eolis.org/gridded-product/), we know that it is generated at\n",
"a monthly temporal resolution.\n",
"\n",
"We can now construct a monthly cumulative elevation change timeseries for glaciers of each surge-type."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def calculate_time_series(gridded_product_data: gpd.GeoDataFrame) -> Tuple[NDArray, NDArray, List[int]]:\n",
" \"\"\"\n",
" Calculate a weighted average elevation change per month and associated errors.\n",
"\n",
" Weights are assigned using the CryoTEMPO-EOLIS gridded product uncertainty.\n",
"\n",
" Parameters\n",
" ----------\n",
" gridded_product_data : gpd.GeoDataFrame\n",
" CryoTEMPO-EOLIS gridded product data.\n",
"\n",
" Returns\n",
" -------\n",
" Tuple[NDArray, NDArray, List[int]]\n",
" Weighted average elevation change elevation change per month, error on weighted average,\n",
" and dates of each timestep.\n",
" \"\"\"\n",
" unique_timestamps = sorted(gridded_product_data.timestamp.unique())\n",
"\n",
" # elevation change is in reference to first gridded product\n",
" reference_dem = gridded_product_data[gridded_product_data.timestamp == unique_timestamps[0]]\n",
"\n",
" # initialise empty lists for time series data\n",
" weighted_mean_timeseries = []\n",
" weighted_mean_errors = []\n",
" for unique_timestamp in unique_timestamps:\n",
" gridded_product_for_timestamp = gridded_product_data[gridded_product_data.timestamp == unique_timestamp]\n",
" merged_gdf = reference_dem.sjoin(df=gridded_product_for_timestamp, how='inner')\n",
" merged_gdf['elevation_difference'] = merged_gdf['elevation_right'] - merged_gdf['elevation_left']\n",
"\n",
" # combine uncertainties for elevation difference\n",
" merged_gdf['elevation_difference_unc'] = np.sqrt(\n",
" merged_gdf['uncertainty_right']**2 + merged_gdf['uncertainty_left']**2)\n",
"\n",
" # calculate average elevation change, weighted by measurement uncertainty\n",
" weighted_mean = (np.sum(merged_gdf['elevation_difference'] / merged_gdf['elevation_difference_unc']**2)\n",
" / np.sum(1 / merged_gdf['elevation_difference_unc']**2))\n",
"\n",
" # calculate weighted average uncertainty\n",
" error = np.sqrt(1 / np.sum(1 / merged_gdf['elevation_difference_unc']**2))\n",
"\n",
" weighted_mean_timeseries.append(weighted_mean)\n",
" weighted_mean_errors.append(error)\n",
"\n",
" dates = [datetime.fromtimestamp(ts) for ts in unique_timestamps]\n",
"\n",
" return np.array(weighted_mean_timeseries), np.array(weighted_mean_errors), dates"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plotting this time series, we see how the glaciers that are known to surge, or possible surging glaciers, are reducing in elevation more than non-surging glaciers."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots(figsize=(10, 5))\n",
"\n",
"for i, (surge_type, gridded_product_data) in enumerate(svalbard_eolis_data.items()):\n",
" # calculate time series for surge type\n",
" weighted_mean_timeseries, weighted_mean_errors, dates = calculate_time_series(gridded_product_data)\n",
" \n",
" # plot time series\n",
" color = cm.CMRmap((i + 1) / (len(surge_type_definitions)))\n",
" ax.plot(dates, weighted_mean_timeseries, lw=1, label=surge_type, c=color)\n",
" ax.fill_between(dates,\n",
" weighted_mean_timeseries - weighted_mean_errors,\n",
" weighted_mean_timeseries + weighted_mean_errors,\n",
" alpha=0.5,\n",
" color=color)\n",
"\n",
"ax.set_xlabel('Date')\n",
"ax.set_ylabel('Cumulative Elevation Change [m]')\n",
"ax.legend(title='Surge type')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Map rate of surface elevation change\n",
"\n",
"Finally, we can calculate the rate of surface elevation change per year by applying a least squares regression to all data points for each pixel.\n",
"\n",
"We can plot this spatially for the different surge-types to see how the glacier elevation is changing across svalbard, and plot the mean rate of elevation change per surge-type."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def calculate_rate_of_change_map(gridded_product_data: gpd.GeoDataFrame) -> gpd.GeoDataFrame:\n",
" \"\"\"\n",
" Apply least squares linear regression to each pixel through time.\n",
"\n",
" Parameters\n",
" ----------\n",
" gridded_product_data : gpd.GeoDataFrame\n",
" CryoTEMPO-EOLIS gridded product data.\n",
"\n",
" Returns\n",
" -------\n",
" gpd.GeoDataFrame\n",
" Rate of surface elevation change per year per pixel.\n",
" \"\"\"\n",
" pixel_groups = gridded_product_data.groupby(by='geometry')\n",
" rate_of_surface_elevation_change_years = []\n",
" for _, pixel_group in pixel_groups:\n",
" pixel_group_df = pixel_group.sort_values(by='timestamp')\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.elevation.to_numpy(),\n",
" rcond=None)[0]\n",
" # m gives the rate of change of elevation in seconds, lets crudely convert it to years for a more meaningful value\n",
" rate_of_surface_elevation_change_years.append(m * 60 * 60 * 24 * 365)\n",
" dhdt_gdf = gpd.GeoDataFrame({'rate_of_surface_elevation_change_years': rate_of_surface_elevation_change_years},\n",
" geometry=[key for key, _ in pixel_groups], crs=gridded_product_data.crs)\n",
" return dhdt_gdf "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, axes = plt.subplots(figsize=(20, 20), nrows=2, ncols=2)\n",
"axes = axes.flatten()\n",
"\n",
"min_x, min_y, max_x, max_y = rgi_data.to_crs(nps_epsg_code).total_bounds\n",
"\n",
"mean_rate_of_change_for_surge_types = []\n",
"for i, (surge_type, gridded_product_data) in enumerate(svalbard_eolis_data.items()):\n",
" # calculate rate of surface elevation change for pixel over time\n",
" dhdt_gdf = calculate_rate_of_change_map(gridded_product_data)\n",
" mean_rate_of_change_for_surge_types.append(dhdt_gdf['rate_of_surface_elevation_change_years'].mean())\n",
"\n",
" # plot rate of surface elevation change on map\n",
" dhdt_gdf.to_crs(epsg=nps_epsg_code).plot(\n",
" ax=axes[i], column='rate_of_surface_elevation_change_years', s=5, cmap='inferno', vmin=-5, vmax=5,\n",
" legend=True, legend_kwds={\"label\": \"Elevation change [m / year]\"})\n",
" \n",
" axes[i].set_xlim(min_x, max_x)\n",
" axes[i].set_ylim(min_y, max_y)\n",
" axes[i].set_xlabel('x [m]')\n",
" axes[i].set_ylabel('y [m]')\n",
" axes[i].set_title(f'Rate of surface elevation change \\n for surge type: {surge_type}')\n",
"\n",
" ctx.add_basemap(axes[i], source=ctx.providers.Esri.WorldImagery, crs=nps_epsg_code, zoom=6)\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, ax = plt.subplots(figsize=(5, 5))\n",
"# round mean rate of change to 2 decimal places\n",
"rects = ax.bar(svalbard_eolis_data.keys(), [round(mean_rate, 2) for mean_rate in mean_rate_of_change_for_surge_types],\n",
" color=[cm.CMRmap((i + 1) / (len(surge_type_definitions))) for i in range(len(svalbard_eolis_data))])\n",
"ax.bar_label(rects, padding=3)\n",
"plt.ylabel('Mean rate of surface elevation change \\n (metres / year)')\n",
"plt.title('Mean Rate of Surface Elevation Change for \\n Different Surge-Type Glaciers in Svalbard')\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
}