shean_ICESat-2_hackweek_tutorial_GeospatialAnalysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Geospatial Data Exploration, Analysis, and Visualization\n",
    "#### A sample of relevant concepts and methods that can be used for ICESat-2 data analysis*\n",
    "*Not intended for billions of ATL03 points\n",
    "\n",
    "ICESat-2 hackweek  \n",
    "June 12, 2020  \n",
    "David Shean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You've used the NSIDC API, pulled your granules, extracted the relevant records and fields, and done some initial QA/QC.  Now what?  \n",
    "\n",
    "Let's work through some Python geospatial data science tools that can help with analysis and visualization of these point data, enabling exploratory data analysis and data-driven discovery.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Objectives\n",
    "\n",
    "* Review fundamental concepts in geospatial analysis (e.g. coordinate systems, projections, datums, geometry types)\n",
    "* Learn basic geospatial data manipulation and exploration with a relatively small, clean ICESat GLAS dataset\n",
    "* Use modern data science tools (pandas) for basic data exploration, interpretation, and visualization\n",
    "* Explore different approaches for spatial aggregation of points - groupby attribute, spatial join by polygon, hex bins\n",
    "* Learn how to sample a raster at point locations\n",
    "* Analyze elevation change using sparse laser altimetry data over mountain glaciers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial prep\n",
    "\n",
    "1. Log on to the event Jupyterhub: https://icesat-2.hackweek.io\n",
    "\n",
    "1. Open a terminal in Jupyterlab (\"+\" icon in upper left corner, then click on the terminal icon)\n",
    "\n",
    "1. Navigate to local directory where you want to store tutorial material (home directory by default):  \n",
    "`cd ~/tutorials`\n",
    "\n",
    "1. If you haven't dont so, clone the repo:  \n",
    "`git clone https://github.com/ICESAT-2HackWeek/geospatial-analysis.git`\n",
    "\n",
    "1. Enter the tutorial directory:  \n",
    "`cd geospatial-analysis`\n",
    "\n",
    "    1. If you cloned earlier in the week, you will need to pull latest changes:\n",
    "`git pull`\n",
    "\n",
    "**Remeber to post all questions to the `#questions` channel on the Slack workspace**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quick Zoom poll:\n",
    "In the Zoom Participants panel, use the Yes (green) or No (red) response to answer the following question: \n",
    "\n",
    "**Have you ever taken a GIS course?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## GIS Basics\n",
    ">A geographic information system (GIS) is a framework for gathering, managing, and analyzing data. Rooted in the science of geography, GIS integrates many types of data. It analyzes spatial location and organizes layers of information into visualizations using maps and 3D scenes. ​With this unique capability, GIS reveals deeper insights into data, such as patterns, relationships, and situations—helping users make smarter decisions.\n",
    "\n",
    "https://www.esri.com/en-us/what-is-gis/overview\n",
    "\n",
    "Primary data types:\n",
    "* Vector - points, lines, polygons; shapefiles\n",
    "    * https://automating-gis-processes.github.io/site/notebooks/L1/geometric-objects.html\n",
    "* Raster - images, gridded data; GeoTiff\n",
    "\n",
    "Concepts (with specific geopandas doc):\n",
    "* Coordinate Systems - map projection and datum: https://geopandas.org/projections.html\n",
    "* Geometry manipulation/operations: \n",
    "    * https://geopandas.org/geometric_manipulations.html\n",
    "    * https://geopandas.org/set_operations.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Scientific Python landscape\n",
    "\n",
    "* Python\n",
    "* Jupyter/iPython\n",
    "* NumPy, Pandas, Matplotlib, SciPy\n",
    "* xarray, scikit-learn\n",
    "\n",
    "One (aging) interpretation of this stack:\n",
    "\n",
    "![2017 Scientific Python Stack](https://devopedia.org/images/article/60/7938.1587985662.jpg)  \n",
    "Slide from Jake VanderPlas’s presentation at PyCon 2017, entitled “The Unexpected Effectiveness of Python in Science.”"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Geospatial Python landscape\n",
    "* GDAL/OGR, GEOS, PROJ\n",
    "* rasterio, fiona, shapely, pyproj\n",
    "* geopandas, cartopy, xarray\n",
    "* rioxarray/salem, PDAL\n",
    "\n",
    "[Insert shiny new figure here]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Many excellent resources available\n",
    "Here are a few...\n",
    "\n",
    "#### Full courses:\n",
    "* https://www.earthdatascience.org/courses/\n",
    "* https://automating-gis-processes.github.io\n",
    "* https://github.com/UW-GDA/gda_course_2020\n",
    "\n",
    "#### Geohackweek: \n",
    "* https://geohackweek.github.io/\n",
    "\n",
    "#### Pangeo:\n",
    "* https://github.com/pangeo-data/pangeo-tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Complementary approaches for ICESat-2 data\n",
    "\n",
    "1. Efficient, scalable processing of huge point datasets\n",
    "    * Mostly NumPy\n",
    "    * Basic array manipulations\n",
    "    * Array size limited by memory - can chunk with Dask\n",
    "2. Higher-level data science - *analysis, interpetation, and visualization*\n",
    "    * NumPy under the hood\n",
    "    * Convenience and flexibility comes with small performance hit\n",
    "    * Labels make life easier (don't need to remember integer indices)\n",
    "\n",
    "Here, we're going to explore #2 with an existing dataset.\n",
    "\n",
    "At the end of the day, most applications just need clean (x,y,z,t) points.\n",
    "\n",
    "As with all things in the *nix/open-source/Python world, there are always multiple approaches that can be used to accomplish the same goals.  The user must decide on an approach based on complexity, time constraints, etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ICESat GLAS Background\n",
    "The NASA Ice Cloud and land Elevation Satellite ([ICESat](https://icesat.gsfc.nasa.gov/icesat/)) was a NASA mission carrying the Geosciences Laser Altimeter System (GLAS) instrument: a space laser, pointed down at the Earth (and unsuspecting Earthlings).  \n",
    "\n",
    "It measured surface elevations by precisely tracking laser pulses emitted from the spacecraft at a rate of 40 Hz (a new pulse every 0.025 seconds).  These pulses traveled through the atmosphere, reflected off the surface, back up through the atmosphere, and into space, where some small fraction of that original energy was received by a telescope on the spacecraft.  The instrument electronics precisely recorded the time when these intrepid photons left the instrument and when they returned.  The position and orientation of the spacecraft was precisely known, so the two-way traveltime (and assumptions about the speed of light and propagation through the atmosphere) allowed for precise forward determination of the spot on the Earth's surface (or cloud tops, as was often the case) where the reflection occurred.  The laser spot size varied during the mission, but was ~70 m in diameter. \n",
    "\n",
    "ICESat collected billions of measurements from 2003 to 2009, and was operating in a \"repeat-track\" mode that sacrificed spatial coverage for more observations along the same ground tracks over time.  One primary science focus involved elevation change over the Earth's ice sheets.  It allowed for early measurements of full Antarctic and Greenland ice sheet elevation change, which offered a detailed look at spatial distribution and rates of mass loss, and total ice sheet contributions to sea level rise.  \n",
    "\n",
    "There were problems with the lasers during the mission, so it operated in short campaigns lasting only a few months to prolong the full mission lifetime.  While the primary measurements focused on the polar regions, many measurements were also collected over lower latitudes, to meet other important science objectives (e.g., estimating biomass in the Earth's forests, observing sea surface height/thickness over time). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 1: Sample GLAS dataset for CONUS\n",
    "A few years ago, I wanted to evaluate ICESat coverage of the Continental United States (CONUS).  The primary application was to extract a set of accurate control points to co-register a large set of high-resolution digital elevation modoels (DEMs) derived from satellite stereo imagery.  I wrote some Python/shell scripts to download, filter, and process all of the GLAH14 granules in parallel ([https://github.com/dshean/icesat_tools](https://github.com/dshean/icesat_tools)).\n",
    "\n",
    "The high-level workflow is here: https://github.com/dshean/icesat_tools/blob/master/glas_proc.py#L24.  These tools processed each HDF5 (H5) file and wrote out csv files containing only the “good” points. These csv files were concatenated to prepare the single input csv (`GLAH14_tllz_conus_lulcfilt_demfilt.csv`) that we will use for this tutorial.  \n",
    "\n",
    "The csv contains ICESat GLAS shots that passed the following filters:\n",
    "* Within some buffer (~110 km) of mapped glacier polygons from the [Randolph Glacier Inventory (RGI)](https://www.glims.org/RGI/)\n",
    "* Returns from exposed bare ground (landcover class 31) or snow/ice (12) according to a 30-m Land-use/Land-cover dataset (2011 NLCD, https://www.mrlc.gov/data?f%5B0%5D=category%3Aland%20cover)\n",
    "* Elevation values within some threshold (200 m) of elevations sampled from an external reference DEM (void-filled 1/3-arcsec [30-m] SRTM-GL1, https://lpdaac.usgs.gov/products/srtmgl1v003/), used to remove spurious points and returns from clouds.\n",
    "* Various other ICESat-specific quality flags (see comments in `glas_proc.py` for details)\n",
    "\n",
    "The final file contains a relatively small subset (~65K) of the total shots in the original GLAH14 data granules from the full mission timeline (2003-2009).  The remaining points should represent returns from the Earth's surface with reasonably high quality, and can be used for subsequent analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wait, I thought this was an ICESat-**2** hackweek?\n",
    "\n",
    "Note that we could (and should!) swap similarly processed/filtered ATL06 points over CONUS for this tutorial.  I did not have time to do this before the hackweek, but it would make a great project (nudge).  After opening this updated file with Pandas, it should mostly be a matter of updating field names throughout the notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Import necessary modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import requests\n",
    "import zipfile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "from shapely.geometry import Point, Polygon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Magic function to enable interactive plotting (zoom/pan) in Jupyter notebook\n",
    "#If running locally, this would be `%matplotlib notebook`, but since we're using Juptyerlab, we use widget\n",
    "%matplotlib widget\n",
    "#%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define path to sample GLAS data\n",
    "glas_fn = 'GLAH14_tllz_conus_lulcfilt_demfilt.csv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Quick check of csv file contents\n",
    "!head $glas_fn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pandas\n",
    "\n",
    "Trust me, you should learn how to use `Pandas`, regardless of your ICESat-2 application.  \n",
    "\n",
    "A significant portion of the Python data science ecosystem is based on Pandas and/or Pandas data models.\n",
    "\n",
    ">pandas is a Python package providing fast, flexible, and expressive data structures designed to make working with \"relational\" or \"labeled\" data both easy and intuitive. It aims to be the fundamental high-level building block for doing practical, real world data analysis in Python. Additionally, it has the broader goal of becoming the most powerful and flexible open source data analysis / manipulation tool available in any language. It is already well on its way towards this goal.\n",
    "\n",
    "https://github.com/pandas-dev/pandas#main-features\n",
    "\n",
    "If you are working with tabular data (rows and columns, like a csv or spreadsheet), especially time series data, please use pandas.\n",
    "* A better way to deal with tabular data, built on top of NumPy arrays\n",
    "* With NumPy, we had to remember which column number (e.g., 3, 4) represented each variable (lat, lon, glas_z, etc)\n",
    "* Pandas allows you to store data with different types, and then reference using more meaningful labels\n",
    "    * NumPy: `glas_np[:,4]`\n",
    "    * Pandas: `glas_df['glas_z']`\n",
    "* A good \"10-minute\" reference with examples: https://pandas.pydata.org/pandas-docs/stable/getting_started/10min.html\n",
    "\n",
    "If you are working with more complex data, like collections of tabular time series data from 100s of met stations or netCDF model output, you can use [`xarray` package](http://xarray.pydata.org/en/stable/), which extends the `pandas` data model to n-dimensions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the csv file with Pandas\n",
    "* Note that pandas has excellent readers for most common file formats: https://pandas.pydata.org/pandas-docs/stable/reference/io.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_df = pd.read_csv(glas_fn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## That was easy.  Let's inspect the returned `DataFrame` object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check data types\n",
    "* Can use the DataFrame `info` method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_df.info()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get the column labels\n",
    "* Can use the DataFrame `columns` attribute"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_df.columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you are new to Python and object-oriented programming, take a moment during the break to consider the [difference between the methods and attributes](https://stackoverflow.com/questions/46312470/difference-between-methods-and-attributes-in-python) of the DataFrame, and how both are accessed. \n",
    "\n",
    "https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.html\n",
    "\n",
    "If this is confusing, ask your neighbor or instructor."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preview records using DataFrame `head` and `tail` methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_df.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute the mean and standard deviation for all values in each column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_df.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_df.std()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Apply a custom function to each column\n",
    "* For this example, let's define a function to compute the Normalized Median Absolute Deviation (NMAD)\n",
    "* https://en.wikipedia.org/wiki/Median_absolute_deviation\n",
    "    * For a normal distribution, this is equivalent to the standard deviation. \n",
    "    * For data containing outliers, it is a more robust representation of variability. \n",
    "* We will then use the Pandas `apply` method to compute the NMAD for all values in each column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def nmad(a, c=1.4826):\n",
    "    return np.median(np.fabs(a - np.median(a))) * c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_df.apply(nmad)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Note: the NMAD function is now distributed with `scipy.stats`: \n",
    "https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.median_absolute_deviation.html*\n",
    "\n",
    "Can also use and `apply` functions imported from modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_df.apply(scipy.stats.median_absolute_deviation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Print quick stats for entire DataFrame with the `describe` method"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_df.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Useful, huh?  Note that the `50%` statistic is the median."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use the Pandas plotting functionality to create a 2D scatterplot of `glas_z` values\n",
    "* https://pandas.pydata.org/pandas-docs/stable/reference/api/pandas.DataFrame.plot.scatter.html\n",
    "* Note that labels and colorbar are automatically plotted!\n",
    "* Can adjust the size of the points using the `s=1` keyword\n",
    "* Can experiment with different color ramps:\n",
    "    * https://matplotlib.org/examples/color/colormaps_reference.html (I prefer `inferno`)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Color ramps\n",
    "Information on how to choose a good colormap for your data: https://matplotlib.org/3.1.0/tutorials/colors/colormaps.html  \n",
    "Another great resource (Thanks @fperez!): https://matplotlib.org/cmocean/  \n",
    "**TL;DR** Don't use `jet`, use a perceptually uniform colormap for linear variables like elevation. Use a diverging color ramp for values where sign is important."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_df.plot(x='lon', y='lat', kind='scatter', c='glas_z', s=1, cmap='inferno');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Experiment by changing the variable represented with the color ramp\n",
    "* Try `decyear` or other columns to quickly visualize spatial distribution of these values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_df.plot(x='lon', y='lat', kind='scatter', c='decyear', s=1, cmap='inferno');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a histogram that shows the number of points vs time (`decyear`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Determine number of bins needed to provide weekly resolution\n",
    "nbins = int(50 * np.ptp(glas_df['decyear'].values))\n",
    "nbins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = glas_df.hist('decyear', bins=nbins)[0,0]\n",
    "ax.set_ylabel('Number of ICESat points')\n",
    "ax.set_title('Weekly usable ICESat point count over CONUS');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we can resolve the distinct campaigns during the ICESat mission (each ~1-2 months long)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a histogram of all `glas_z` elevation values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Determine number of bins needed for 20 m elevation bins\n",
    "binwidth = 20\n",
    "nbins = int(np.ptp(glas_df['glas_z'].values)/binwidth)\n",
    "nbins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_df.hist('glas_z', bins=nbins);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wait a minute...negative elevations!?  Who calibrated this thing? C'mon NASA."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A note on vertical datums\n",
    "\n",
    "Note that some elevations are less than 0 m.  How can this be?\n",
    "\n",
    "The `glas_z` values are height above (or below) the WGS84 ellipsoid. This is not the same vertical datum as mean sea level (roughly approximated by a geoid model).\n",
    "\n",
    "A good resource explaining the details: https://vdatum.noaa.gov/docs/datums.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How many GLAS points have a negative `glas_z` value?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_df['glas_z'] < 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(glas_df['glas_z'] < 0).value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check spatial distribution of points below 0 (height above WGS84 ellipsoid)\n",
    "* Create a scatterplot only using points with negative values\n",
    "* Adjust the color ramp bounds to bring out more detail for these points\n",
    "    * hint: see the `vmin` and `vmax` arguments for the `plot` function\n",
    "* What do you notice about these points? (may be tough without more context, like coastlines and state boundaries or a tiled basemap - we'll learn how to incorporate these later)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_df[glas_df['glas_z'] < 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_df[glas_df['glas_z'] < 0].plot(x='lon', y='lat', kind='scatter', c='dem_z', s=1, cmap='inferno', vmin=-30, vmax=0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Geoid offset grids\n",
    "PROJ and GDAL use raster grids of the horizontal and vertical offsets between the WGS84 ellipsoid and different geoid models.\n",
    "\n",
    "These grids are now hosted on the cloud: https://cdn.proj.org/, which enables on-the-fly transformations with GDAL!\n",
    "\n",
    "Here's the EGM96 offset grid.  It shows the height differene between the WGS84 ellipsoid (simple shape model of the Earth) and the EGM96 geoid, which approximates a geopotential (gravitational) surface, approximately euivalent to mean sea level \n",
    "\n",
    "![EGM96 geoid offset grid](https://raw.githubusercontent.com/UW-GDA/gda_course_2020/master/resources/sample_img/egm96_offset.png)  \n",
    "(Source GeoTiff: https://cdn.proj.org/us_nga_egm96_15.tif)\n",
    "\n",
    "\n",
    "Note values for CONUS.  \n",
    "\n",
    "A lot of the points with elevation < 0 m in the above plot are near coastal sites, roughly near mean sea level.  We see that the geoid offset (difference between WGS84 ellipsoid and EGM96 geoid in this case) for CONUS is roughly -20 m.  So the ICESat GLAS point elevations near the coast are roughly -20 m relative to the ellipsoid, even though they are 0 m relative to the geoid (approximately mean sea level).  These concepts around datums can be a bit confusing, so please ask questions.\n",
    "\n",
    "Note that ICESat-2 includes elevation values relative to the EGM2008 geoid model to provide orthometric heights - a more recent, more accurate, more detailed geoid model compared to EGM96."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 2: More Pandas, Outlier Removal, Groupby"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute the elevation difference between ICESat `glas_z` and SRTM `dem_z` values\n",
    "\n",
    "Earlier, I mentioned that I had sampled the SRTM DEM for each GLAS point.  See Appendix for an example of how to do this with rasterio.  \n",
    "\n",
    "For now, let's use the previously sampled values, compute an elevation difference and store in a new column in our DataFrame called `glas_srtm_dh`\n",
    "\n",
    "Remember the order of this calculation (if the difference values are negative, which dataset is higher elevation?)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_df['glas_srtm_dh'] = glas_df['glas_z'] - glas_df['dem_z']\n",
    "glas_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute the time difference between ICESat point timestamp and the SRTM timestamp\n",
    "* Store in a new column named `glas_srtm_dt`\n",
    "* The SRTM data were collected between February 11-22, 2000\n",
    "    * Can assume a constant decimal year value of 2000.112 for now"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#February 11-22, 2000\n",
    "srtm_decyear = 2000.112\n",
    "glas_df['glas_srtm_dt'] = glas_df['decyear'] - srtm_decyear\n",
    "glas_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute *apparent* annualized elevation change rate ($\\frac{dh}{dt}$ in meters per year) from these new columns\n",
    "* This will be rate of change between the SRTM timestamp (2000) and each GLAS point timestamp (2003-2009)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_df['glas_srtm_dhdt'] = glas_df['glas_srtm_dh']/glas_df['glas_srtm_dt']\n",
    "glas_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a scatterplot of the difference values\n",
    "* Use a symmetrical `RdBu` (Red to Blue) color ramp\n",
    "* Set the color ramp limits using `vmin` and `vmax` keyword arguments to be symmetrical about 0 z"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = glas_df.plot(x='lon', y='lat', kind='scatter', c='glas_srtm_dh', s=1, cmap='RdBu', vmin=-10, vmax=10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = glas_df.plot(x='lon', y='lat', kind='scatter', c='glas_srtm_dh', s=1, cmap='RdBu', vmin=-50, vmax=50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute some descriptive statistics\n",
    "* Why might we have a non-zero mean/median difference?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(glas_df['glas_srtm_dh'].mean())\n",
    "print(glas_df['glas_srtm_dh'].std())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(glas_df['glas_srtm_dh'].median())\n",
    "print(nmad(glas_df['glas_srtm_dh']))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a histogram of the difference values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots()\n",
    "glas_df.hist('glas_srtm_dh', ax=ax, bins=128, range=(-10,10));\n",
    "ax.axvline(0, color='k')\n",
    "ax.axvline(glas_df['glas_srtm_dh'].median(), color='r');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a scatterplot of elevation difference `glas_srtm_dh` values vs elevation values\n",
    "* `glas_srtm_dh` should be on the y-axis\n",
    "* `glas_z` values on the x-axis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = glas_df.plot('glas_z', 'glas_srtm_dh', kind='scatter', s=1)\n",
    "#Add a horizontal line at 0\n",
    "ax.axhline(0, color='k', lw=0.5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Remove outliers\n",
    "The initial filter in `glas_proc.py` removed GLAS points with absolute elevation difference >200 m compared to the SRTM elevations.  We expect most real elevation change signals to be less than this for the given time period.  But clearly some outliers remain.\n",
    "\n",
    "Let's design and apply a simple filter that removes outliers.  One option is to define outliers as values outside some absolute threshold. Can set this threshold as some multiple of the standard deviation (e.g., `3*std`). Can also use quantile or percentile values for this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Mean difference:\", glas_df['glas_srtm_dh'].mean())\n",
    "thresh = 3.5 * glas_df['glas_srtm_dh'].std()\n",
    "print(\"3.5 * std:\", thresh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = (glas_df['glas_srtm_dh'] - glas_df['glas_srtm_dh'].mean()).abs() <= thresh\n",
    "glas_df_fltr = glas_df[idx]\n",
    "print(\"Number of points before filter:\", glas_df.shape[0])\n",
    "print(\"Number of points after filter:\", glas_df_fltr.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clim = thresh\n",
    "f, axa = plt.subplots(1,2, figsize=(10,4))\n",
    "#Outliers plotted in black\n",
    "glas_df.plot(ax=axa[1], x='glas_z', y='glas_srtm_dh', kind='scatter', s=1, color='k', label='Outliers')\n",
    "glas_df_fltr.plot(ax=axa[0], x='lon', y='lat', kind='scatter', c='glas_srtm_dh', s=1, cmap='RdBu', vmin=-clim, vmax=clim)\n",
    "glas_df_fltr.plot(ax=axa[1], x='glas_z', y='glas_srtm_dh', kind='scatter', s=1, c='orange', label='Inliers')\n",
    "glas_df[~idx].plot(ax=axa[0], x='lon', y='lat', kind='scatter', color='k', s=1, legend=False)\n",
    "axa[1].axhline(0,color='k')\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Active remote sensing sanity check\n",
    "\n",
    "Even after removing outliers, there are still some big differences between the SRTM and GLAS elevation values.  \n",
    "\n",
    "* Do you see systematic differences between the glas_z and dem_z values?\n",
    "* Any clues from the scatterplot? (e.g., do some tracks (north-south lines of points) display systematic bias?)\n",
    "* Brainstorm some ideas about what might be going on here.  Think about the nature of each sensor:\n",
    "    * ICESat was a Near-IR laser (1064 nm wavelength) with a big ground spot size (~70 m in diameter)\n",
    "        * Timestamps span different seasons between 2003-2009\n",
    "    * SRTM was a C-band radar (5.3 GHz, 5.6 cm wavelength) with approximately 30 m ground sample distance (pixel size)\n",
    "        * Timestamp was February 2000\n",
    "        * Data gaps (e.g., radar shadows, steep slopes) were filled with ASTER GDEM2 composite, which blends DEMs acquired over many years ~2000-2014\n",
    "* Consider different surfaces and how the laser/radar footprint might be affected:\n",
    "    * Flat bedrock surface\n",
    "    * Dry sand dunes\n",
    "    * Steep montain topography like the Front Range in Colorado  \n",
    "    * Dense vegetation of the Hoh Rainforest in Olympic National Park"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Let's check to see if differences are due to our land-use/land-cover (LULC) classes\n",
    "* Determine the unique values in the `lulc` column (hint: see the `value_counts` method)\n",
    "* In the introduction, I mentioned that I initially preserved only two classes for these GLAS points (12 - snow/ice, 31 - barren land), so this isn't going to help us over forests:\n",
    "    * https://www.mrlc.gov/data/legends/national-land-cover-database-2011-nlcd2011-legend"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_df['lulc'].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use Pandas `groupby` to compute stats for the LULC classes\n",
    "* https://pandas.pydata.org/pandas-docs/stable/user_guide/groupby.html\n",
    "* This is one of the most powerful features in Pandas, efficient grouping and analysis based on attributes\n",
    "* Compute mean, median and std of the difference values (glas_z - dem_z) for each LULC class\n",
    "* Do you see a difference between values over glaciers vs bare rock?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_df[['glas_srtm_dh', 'lulc']].groupby('lulc').mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Statistics for full dataset\n",
    "glas_df[['glas_srtm_dh', 'lulc']].groupby('lulc').agg(['mean', 'std', 'median', nmad])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Statistics for filtered dataframe (after removing outliers)\n",
    "glas_df_fltr[['glas_srtm_dh', 'lulc']].groupby('lulc').agg(['mean', 'std', 'median', nmad])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for title, group in glas_df_fltr.groupby('lulc'):\n",
    "    group.plot(x='lon', y='lat', kind='scatter', c='glas_srtm_dh', s=1, cmap='RdBu', vmin=-10, vmax=10, title='LULC %i' % title)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 3: Geopandas Intro\n",
    "\n",
    "pandas is great, but what if we want to do some geospatial operations - like reproject our points or compute the intersection between Point and Polygon features?\n",
    "\n",
    "Enter Geopandas - all the great things about pandas, plus geo! (http://geopandas.org/).\n",
    "\n",
    ">\"GeoPandas is an open source project to make working with geospatial data in python easier. GeoPandas extends the datatypes used by pandas to allow spatial operations on geometric types. Geometric operations are performed by shapely. Geopandas further depends on fiona for file access and descartes and matplotlib for plotting.\"\n",
    "\n",
    ">\"GeoPandas enables you to easily do operations in python that would otherwise require a spatial database such as PostGIS.\"\n",
    "\n",
    "Under the hood, GeoPandas is `pandas` plus some other core geospatial packages:\n",
    "* `shapely` for geometry operations (https://shapely.readthedocs.io/en/stable/manual.html)\n",
    "* `fiona` for reading/writing GIS file formats (https://fiona.readthedocs.io/en/latest/manual.html)\n",
    "* `pyproj` for projections and coordinate system transformations (http://pyproj4.github.io/pyproj/stable/)\n",
    "\n",
    "Under those hoods are lower-level geospatial libraries (GEOS, GDAL/OGR, PROJ4) that provide a foundation for most GIS software (open-source and commercial).  I encourage you to explore these - I guarantee you will learn something valuable.\n",
    "\n",
    "For now, let's explore some basic geopandas functionality."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convert pandas `DataFrame` to geopandas `GeoDataFrame`\n",
    "To do this, we need to create a new column containing standardized `geometry` objects (e.g., `Point`, `Polygon`) for each record in the DataFrame. \n",
    "\n",
    "https://geopandas.readthedocs.io/en/latest/gallery/create_geopandas_from_pandas.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Convert the Pandas `DataFrame` to a GeoPandas `GeoDataFrame`\n",
    "* https://geopandas.readthedocs.io/en/latest/gallery/create_geopandas_from_pandas.html\n",
    "* Careful about lon and lat order!\n",
    "* Define coordinate reference system (4326 is geographic lat/lon on WGS84 Ellispoid) \n",
    "    * https://www.earthdatascience.org/courses/earth-analytics/spatial-data-r/intro-to-coordinate-reference-systems/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf = gpd.GeoDataFrame(glas_df, geometry=gpd.points_from_xy(glas_df['lon'], glas_df['lat']), crs='EPSG:4326')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(glas_gdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(glas_df)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Same `head()` method works on a GeoDataFrame! inherited from Pandas DataFrame\n",
    "glas_gdf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#But we have additional attributes/methods in a GeoDataFrame\n",
    "glas_gdf.crs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a quick 2D scatterplot\n",
    "\n",
    "Like a Pandas DataFrame, a GeoDataFrame has convenience plotting function that is built on matlplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, looks like a scatterplot.  But let's plot the elevation values with a color ramp.  \n",
    "To do this, just specify the column name as the first argument to `plot`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf.plot('glas_z', markersize=1, cmap='inferno', legend=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load and plot state polygons for context\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hmmm, let's see.  Two choices:\n",
    "1. We could go to ESRI or the U.S. Census website, identify and download a shapefile, unzip 4+ files, copy/paste the appropriate \\*.shp filename into the notebook.  Wait, how can I download on a remote server?  OK, maybe run something like `wget http://...`, unzip, provide absolute path  \n",
    "*- OR -*\n",
    "2. Give geopandas a url that points to a GeoJSON file somewhere on the web, and read dynamically\n",
    "\n",
    "Yeah, let's go with #2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use the US States 5M GeoJSON here: http://eric.clst.org/tech/usgeojson/\n",
    "\n",
    "We've heard GeoJSON mentioned a few times this week.  It's a great format.  If you are unfamiliar: https://en.wikipedia.org/wiki/GeoJSON"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#1:5000000 scale polygons\n",
    "states_url = 'http://eric.clst.org/assets/wiki/uploads/Stuff/gz_2010_us_040_00_5m.json'\n",
    "#1:500000 scale polygons (larger file, more vertices)\n",
    "#states_url = 'http://eric.clst.org/assets/wiki/uploads/Stuff/gz_2010_us_040_00_500k.json'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "states_gdf = gpd.read_file(states_url)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inspect State GeoDataFrame\n",
    "Note that some geometry entries are tuples of POLYGON objects - these are states with islands or rings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "states_gdf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "states_gdf.crs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Limit to Lower48"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = states_gdf['NAME'].isin(['Alaska','Puerto Rico','Hawaii'])\n",
    "states_gdf = states_gdf[~idx]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = glas_gdf.plot('glas_z', markersize=1, cmap='inferno', legend=True)\n",
    "states_gdf.plot(ax=ax, facecolor='none', edgecolor='0.5', linewidth=0.5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RGI glacier polygons\n",
    "\n",
    "Let's grab some glacier outline poygons from the Randolph Glacier Inventory (RGI) v6.0: https://www.glims.org/RGI/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quick Zoom poll:\n",
    "In the Zoom Participants panel, use the Yes (green) or No (red) response to answer the following question: \n",
    "\n",
    "**Have you ever worked with RGI polygons for glacier and/or ice caps?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Fetch the zip file for Region 02 (Western North America)\n",
    "rgi_zip_fn = '02_rgi60_WesternCanadaUS.zip'\n",
    "\n",
    "if not os.path.exists(rgi_zip_fn):\n",
    "    url = 'https://www.glims.org/RGI/rgi60_files/' + rgi_zip_fn\n",
    "    myfile = requests.get(url)\n",
    "    open(rgi_zip_fn, 'wb').write(myfile.content)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Unzip the file\n",
    "with zipfile.ZipFile(rgi_zip_fn, 'r') as zip_ref:\n",
    "    zip_ref.extractall('rgi')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Specify the shapefile filename\n",
    "rgi_fn = 'rgi/02_rgi60_WesternCanadaUS.shp'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load RGI shapefile using Geopandas\n",
    "* Very easy with `read_file()` meethod"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_gdf = gpd.read_file(rgi_fn)\n",
    "rgi_gdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#By default a new integer index is created.  Can use the RGI ID as our index\n",
    "#rgi_gdf = rgi_gdf.set_index('RGIId')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a quick plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_gdf.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Update our quick plot of RGI polygons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))\n",
    "ax = world.plot(alpha=0.2)\n",
    "states_gdf.plot(ax=ax, facecolor='none', edgecolor='0.5', linewidth=0.5)\n",
    "rgi_gdf.plot(ax=ax, edgecolor='k', linewidth=0.5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clip RGI polygons to CONUS\n",
    "GeoPandas makes spatial selection easy.  \n",
    "\n",
    "We'll have two options: 1) using a bounding box, and 2) using an arbitrary polygon."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Bounding box"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf.total_bounds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Define min/max variables for each dimension"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xmin, ymin, xmax, ymax = glas_gdf.total_bounds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Create new GeoDataFrame from output of simple spatial filter with GeoPandas `cx` function\n",
    "* https://geopandas.org/indexing.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_gdf_conus = rgi_gdf.cx[xmin:xmax, ymin:ymax]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Quick plot to verify"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_gdf_conus.plot(edgecolor='k', linewidth=0.5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. Clip points to arbitrary Polygon geometry\n",
    "\n",
    "Let's define a Polygon around our area of interest (the GLAS points) \n",
    "\n",
    "To do this, we'll first take the union of our ~65K GLAS points, and then compute the convex hull\n",
    "\n",
    "This will return a Polygon geometry object, which renders nicely in the Jupyter notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_chull = glas_gdf.unary_union.convex_hull"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Check the type\n",
    "type(glas_gdf_chull)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Preview geometry\n",
    "Note that geometry objects (points, lines, polygons, etc.) will render directly in the Jupyter notebook!  Great for quick previews."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_chull"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(glas_gdf_chull)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Compute intersection between all RGI polygons and the convex hull\n",
    "Use the GeoDataFrame `intersects()` function.  \n",
    "This will return a Boolean DataSeries, True if points intersect the polygon, False if they do not"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_gdf_idx = rgi_gdf.intersects(glas_gdf_chull)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_gdf_idx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Extract records with True for the intersection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Number of RGI polygons before:\",rgi_gdf.shape[0])\n",
    "rgi_gdf_conus = rgi_gdf[rgi_gdf_idx]\n",
    "print(\"Number of RGI polygons after:\", rgi_gdf_conus.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Quick plot to verify\n",
    "Note latitude range"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_gdf_conus.plot(edgecolor='k', linewidth=0.5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 4: Reprojection and Coordinate Systems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "All of the above examples used standard geodetic lat/lon coordinate system (EPSG:4326).  This is fine for global analyses and basic visualization.  But remember that the width of a degree of longitude varies with latitude (~111 km near equator, ~0 m near pole), so our plots have scaling and aspect ratio issues.\n",
    "\n",
    "We need to choose a map projection that is appropriate for our data. This decision is important for visualization, but is also critical for precise quantitative analysis. For example, if you want to compute area or volume change, you should use an equal-area projection. If you want to calculate distances between two points, you should use an equidistant projection.\n",
    "\n",
    "https://www.axismaps.com/guide/general/map-projections/\n",
    "\n",
    "Sadly, there is no \"perfect\" projection. You, as the mapmaker or data analyst, are responsible for choosing a projection with the right characteristics for your purposes. Let's explore a bit further, and we'll revisit some general guidelines later."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use GeoPandas to reproject your GeoDataFrame\n",
    "* Use the very convenient `to_crs()` method to reproject: https://geopandas.org/projections.html\n",
    "* Start by reprojecting the points to a Universal Transverse Mercator (UTM), Zone 11N (EPSG:32611)\n",
    "* Store the output as a new GeoDataFrame\n",
    "* Do a quick `head()` and note the new values in the `geometry` column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "utm_crs = 'EPSG:32611'\n",
    "glas_gdf_utm = glas_gdf.to_crs(utm_crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_utm.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_utm.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a new plot of the reprojected points\n",
    "* Note the new coordinate system origin (0,0), units, and aspect ratio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = glas_gdf_utm.plot(column='glas_z', cmap='inferno', markersize=1, legend=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Excellent, but what did we just do?\n",
    "\n",
    "Under the hood, GeoPandas used the `pyproj` library (a Python API for PROJ) to transform each point from one coordinate system to another coordinate system.  \n",
    "\n",
    "I guarantee that you've all done this kind of thing before, you may just not remember it or recognize it in this context. See: https://en.wikipedia.org/wiki/List_of_common_coordinate_transformations\n",
    "\n",
    "In 2D, transforming (x,y) coordinates between different projections (e.g., WGS84 vs. UTM) on the same reference ellipsoid is pretty straightforward.  Things start to get more complicated when you include different ellipsoid models, horizontal/vertical datums, etc.  Oh, also the Earth's surface is not static - plate tectonics make everything more complicated, as time becomes important, and transformations must include a \"kinematic\" component.  \n",
    "\n",
    "Fortunately, the `PROJ` library (https://proj.org/about.html) has generalized much of the complicated math for geodetic coordinate transformations.  It's been under development since the 1980s, and our understanding of the Earth's shape and plate motions has changed dramatically in that time period.  So, still pretty complicated, and there are different levels of support/sophistication in the tools/libraries that use `PROJ`.\n",
    "\n",
    "We aren't going to get into the details here, but feel free to take a look at the Transformations section here to get a sense of how this is actually accomplished: https://proj.org/operations/index.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Let's define a custom projection for the Western U.S.\n",
    "\n",
    "The UTM projection we used above is not the best choice for our data, which actually span 4 UTM zones:\n",
    "https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system#/media/File:Utm-zones-USA.svg. \n",
    "\n",
    "We used Zone 11N, which means that our map will have limited distortion within that zone centered on -110°W.  But distortion will increase beyond the width of the -114° to -108°W zone definition.\n",
    "\n",
    "Let's instead use a custom Albers Equal Area (AEA) projection to minimize distoration over the full spatial extent of our GLAS points\n",
    "\n",
    "To do this, we'll define a PROJ string (https://proj.org/usage/quickstart.html?highlight=definition), which can be interpreted by most Python geopackages (like `pyproj`).\n",
    "\n",
    "The Albers Equal Area projection requires two standard parallels: https://proj.org/operations/projections/aea.html.  Here, we will also specify the center latitude and center longitude for the coordinate system origin."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Define a custom Albers Equal Area proj string `'+proj=aea...'`\n",
    "    * https://en.wikipedia.org/wiki/Albers_projection\n",
    "    * PROJ reference, with example: https://proj.org/operations/projections/aea.html\n",
    "    * Here is another sample proj string example for an Albers Equal Area projection of the U.S. (note that this uses GRS80 ellipsoid and NAD83 datum): http://spatialreference.org/ref/esri/102003/proj4/\n",
    "        * For your string, please use WGS84 ellipsoid (see the proj doc for aea in the above link)\n",
    "* Use the center longitude and center latitude you calculated earlier\n",
    "* Define the two standard parallels (lines of latitude) based on the range of your points\n",
    "    * Scale is true along these parallels, and distortion increases as you move away from these two parallels\n",
    "    * One approach would be to use min and max latitude from the `total_bounds` extent computed earlier\n",
    "        * This is fine, but note that this could lead to additional distortion near your center latitude\n",
    "        * Extra Credit: figure out how to place them slightly inside your min and max latitude to minimize distortion across the entire latitude range\n",
    "* Use Python string formatting to dynamically create your proj string (don't just hardcode your values, but substitute variables in the string)\n",
    "* Print the final proj string"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get bounding box (extent) and center (lon, lat) of GLAS points\n",
    "* See GeoPandas API reference. In this case, you want the total_bounds attribute: http://geopandas.org/reference.html#geopandas.GeoSeries.total_bounds\n",
    "* Center can be calculated from the min/max extent values in each dimension"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_extent = glas_gdf.total_bounds\n",
    "glas_extent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_center = glas_extent.reshape((2,2)).mean(axis=0)\n",
    "glas_center"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create the PROJ string"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clon = glas_center[0]\n",
    "clat = glas_center[1]\n",
    "#Standard parallels at min/max values\n",
    "#p1 = glas_extent[1]\n",
    "#p2 = glas_extent[3]\n",
    "#Standard parallels within the point extent\n",
    "dp = (glas_extent[3] - glas_extent[1])*0.67\n",
    "p1 = clat - dp/2.\n",
    "p2 = clat + dp/2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "proj_str_aea = '+proj=aea +lat_1={:.2f} +lat_2={:.2f} +lat_0={:.2f} +lon_0={:.2f}'.format(p1, p2, clat, clon)\n",
    "proj_str_aea"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reproject the GLAS points to the custom projection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_aea = glas_gdf.to_crs(proj_str_aea)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_aea.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bbox = glas_gdf_aea.total_bounds\n",
    "bbox"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create scatter plots for each of the three projections (WGS84, UTM, and custom AEA)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, axa = plt.subplots(1,3, figsize=(10,2.5))\n",
    "glas_gdf.plot(ax=axa[0], column='glas_z', cmap='inferno', markersize=1, legend=True)\n",
    "axa[0].set_title('WGS84')\n",
    "axa[0].grid()\n",
    "glas_gdf_utm.plot(ax=axa[1], column='glas_z', cmap='inferno', markersize=1, legend=True)\n",
    "axa[1].set_title(utm_crs)\n",
    "axa[1].grid()\n",
    "glas_gdf_aea.plot(ax=axa[2], column='glas_z', cmap='inferno', markersize=1, legend=True)\n",
    "axa[2].set_title('Albers Equal-area')\n",
    "axa[2].grid()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Umm, they kind of look the same.  Why are we wasting time on this?\n",
    "\n",
    "* Note the location of the origin for each coordinate system\n",
    "    * The (0,0) should be near the center of your points for the AEA projection\n",
    "    * Where's the origin for the UTM projection?\n",
    "* Note how azimuth and distances are distorted around edges of the plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 5: Spatial Aggregation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hexbin plots\n",
    "Hexbin plots are nice option to visualize the spatial distribution of binned point density or other metric (e.g., median elevation) on a regular grid.  Later in the week, we'll learn about other interpolation and gridding options.\n",
    "\n",
    "Hexagons are preferable over a standard square/rectangular grid: https://pro.arcgis.com/en/pro-app/tool-reference/spatial-statistics/h-whyhexagons.htm\n",
    "\n",
    "Also, see Sarah Battersby's publication on the topic (and appreciate one of the better journal article titles I've seen): https://doi.org/10.1080/15230406.2016.1180263\n",
    "\n",
    "Here are some resources on generating hexbins using Python and matplotlib:\n",
    "* https://matplotlib.org/api/_as_gen/matplotlib.pyplot.hexbin.html\n",
    "* http://darribas.org/gds15/content/labs/lab_09.html\n",
    "\n",
    "Note: an equal-area projection is also a good idea for a hexbin plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nbins = 40"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#To calculate number of bins dynamically with specified bin spacing in meters\n",
    "#bin_width = 27000 #meters\n",
    "#nbins_x = int(np.ceil(np.abs(bbox[2] - bbox[0])/bin_width))\n",
    "#nbins_y = int(np.ceil(np.abs(bbox[3] - bbox[1])/bin_width))\n",
    "#print(nbins_x, nbins_y)\n",
    "#nbins = nbins_x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f,ax = plt.subplots(figsize=(8,6))\n",
    "hb = ax.hexbin(glas_gdf_aea.geometry.x, glas_gdf_aea.geometry.y, gridsize=nbins, cmap='inferno', bins='log', alpha=0.6, mincnt=1)\n",
    "plt.colorbar(hb, label='Point Count')\n",
    "states_gdf.to_crs(proj_str_aea).plot(ax=ax, facecolor='none', edgecolor='black');\n",
    "#Limit plot extent to GLAS point bounds (note: ax.autoscale(False) doesn't work here)\n",
    "ax.set_xlim(bbox[[0,2]]);\n",
    "ax.set_ylim(bbox[[1,3]]);\n",
    "ax.set_title('GLAS Bin Point Count');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f,ax = plt.subplots(figsize=(8,6))\n",
    "hb = ax.hexbin(glas_gdf_aea.geometry.x, glas_gdf_aea.geometry.y, C=glas_gdf_aea['glas_z'], \\\n",
    "               reduce_C_function=np.median, gridsize=nbins, cmap='inferno', alpha=0.6)\n",
    "plt.colorbar(hb, ax=ax, label='Elevation (m HAE)')\n",
    "states_gdf.to_crs(proj_str_aea).plot(ax=ax, facecolor='none', edgecolor='black');\n",
    "ax.set_xlim(bbox[[0,2]]);\n",
    "ax.set_ylim(bbox[[1,3]]);\n",
    "ax.set_title('GLAS Bin Median Elevation');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f,ax = plt.subplots(figsize=(8,6))\n",
    "hb = ax.hexbin(glas_gdf_aea.geometry.x, glas_gdf_aea.geometry.y, C=glas_gdf_aea['glas_srtm_dhdt'], \\\n",
    "               reduce_C_function=np.median, gridsize=nbins, cmap='RdBu', alpha=0.6, vmin=-3, vmax=3)\n",
    "plt.colorbar(hb, ax=ax, label='Elevation Change (m/yr)')\n",
    "states_gdf.to_crs(proj_str_aea).plot(ax=ax, facecolor='none', edgecolor='black');\n",
    "ax.set_xlim(bbox[[0,2]]);\n",
    "ax.set_ylim(bbox[[1,3]]);\n",
    "ax.set_title('Median Elevation Change (m/yr), GLAS - SRTM');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Merge GLAS points with RGI polygons\n",
    "Let's see if we can answer the following question:\n",
    "\n",
    "*Can we identify CONUS glacier surface elevation change that occurred between SRTM (2000) and GLAS (2003-2009) data collection?*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Earlier we computed some statistics for the full CONUS GLAS sample and hex bins.  Now let's analyze the GLAS points that intersect each RGI glacier polygon. \n",
    "\n",
    "One approach would be to loop through each glacier polygon, and do an intersection operation with all points.  But this is inefficient, and doesn't scale well.  It is much more efficient to do a spatial join between the points and the polygons, then groupby and aggregate to compute the relevant statistics for all points that intersect each glacier polygon.\n",
    "\n",
    "You may have learned how to perform a join or spatial join in a GIS course.  So, do we need to open ArcMap or QGIS here?  Do we need a full-fledged spatial database like PostGIS?  No!  GeoPandas has you covered.\n",
    "\n",
    "* Start by reviewing the Spatial Join documentation here: http://geopandas.org/mergingdata.html\n",
    "* Use the geopandas `sjoin` method: http://geopandas.org/reference/geopandas.sjoin.html"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## First, we need to make sure all inputs have the same projection\n",
    "* Reproject the RGI polygons to match our point CRS (custom Albers Equal-area)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_gdf_conus_aea = rgi_gdf_conus.to_crs(glas_gdf_aea.crs)\n",
    "states_gdf_aea = states_gdf.to_crs(glas_gdf_aea.crs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optional: isolate relevant columns to simplify our output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_aea.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_gdf_conus_aea.columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_col = ['glas_z', 'glas_srtm_dhdt', 'geometry']\n",
    "rgi_col = ['RGIId', 'Area', 'Name', 'geometry']\n",
    "\n",
    "glas_gdf_aea_lite = glas_gdf_aea[glas_col]\n",
    "rgi_gdf_conus_aea_lite = rgi_gdf_conus_aea[rgi_col]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_aea_lite"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now try a spatial join between these two \n",
    "* Use the GLAS points as the \"left\" GeoDataFrame and the RGI polygons as the \"right\" GeoDataFrame\n",
    "* Start by using default options (`op='intersects', how='inner'`)\n",
    "* Note the output geometry type and columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_aea_rgi = gpd.sjoin(glas_gdf_aea_lite, rgi_gdf_conus_aea_lite)\n",
    "glas_gdf_aea_rgi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check number of records"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Number of RGI polygons before:\", rgi_gdf_conus_aea.shape[0])\n",
    "print(\"Number of GLAS points before:\", glas_gdf_aea.shape[0])\n",
    "print(\"Number of GLAS points that intersect RGI polygons:\", glas_gdf_aea_rgi.shape[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check number of GLAS points per RGI polygon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_aea_rgi['RGIId'].value_counts()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Which glacier has the greatest number of points?\n",
    "\n",
    "Some notes on indexing and selecting from Pandas DataFrame: http://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html\n",
    "\n",
    "https://www.shanelynn.ie/select-pandas-dataframe-rows-and-columns-using-iloc-loc-and-ix/\n",
    "\n",
    "Here we'll use the `iloc` function to pull out the record for the RGIId key with the highest point count."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "label = glas_gdf_aea_rgi['RGIId'].value_counts().index[0]\n",
    "print(label)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_maxcount = rgi_gdf_conus_aea[rgi_gdf_conus_aea['RGIId'] == label].iloc[0]\n",
    "rgi_maxcount"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_maxcount.geometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ipyleaflet plot\n",
    "\n",
    "OK, great, but where is this glacier?  Let's plot on an interactive ipyleaflet map.\n",
    "\n",
    "Note that leaflet uses tiled basemaps: https://en.wikipedia.org/wiki/Tiled_web_map \n",
    "\n",
    "Default projection is Web Mercator (EPSG:3857): https://en.wikipedia.org/wiki/Web_Mercator_projection.  This works well for lower latitudes, but not the polar regions.\n",
    "\n",
    "FYI, `folium` provides similar functionality outside of the iPython/Jupyter stack: https://python-visualization.github.io/folium/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipyleaflet import Map, Marker, basemaps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Look at all of the basemap options!\n",
    "basemaps.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "center = (rgi_maxcount['CenLat'], rgi_maxcount['CenLon'])\n",
    "basemap = basemaps.Stamen.Terrain\n",
    "m = Map(center=center, zoom=12, basemap=basemap)\n",
    "#label=rgi_gdf_conus_aea.loc[label]['Name']\n",
    "marker = Marker(location=center, draggable=True)\n",
    "m.add_layer(marker);\n",
    "display(m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot points and RGI polygons, Zoom in on WA state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wa_state = states_gdf_aea.loc[states_gdf_aea['NAME'] == 'Washington']\n",
    "#wa_geom = wa_state.iloc[0].geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wa_bbox = wa_state.total_bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clim = (-3.0, 3.0)\n",
    "ax=states_gdf_aea.plot(facecolor='none', edgecolor='black', linewidth=0.5);\n",
    "rgi_gdf_conus_aea.plot(ax=ax, edgecolor='k', lw=0.5, alpha=0.1);\n",
    "glas_gdf_aea.plot(ax=ax, column='glas_srtm_dhdt', cmap='RdBu', markersize=1, legend=True, vmin=clim[0], vmax=clim[1])\n",
    "ax.set_xlim(wa_bbox[[0,2]]);\n",
    "ax.set_ylim(wa_bbox[[1,3]]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Groupby and Aggregate\n",
    "\n",
    "OK, so we know that our sampling isn't great and our dh/dt values are noisy.  But we're here to learn some core concepts and tools, so let's compute some statistics for each glacier anyway.  Hopefully you'll see the value of these operations, and be able to reproduce in the future.\n",
    "\n",
    "We can use the Pandas Groupby functionality to group GLAS points for each RGI polygon, and then aggregate using different functions (e.g., mean, std) for different attributes (e.g., 'glas_z', 'glas_srtm_dhdt').\n",
    "\n",
    "This concept can feel a bit abstract at first, but it is very powerful.  \n",
    "\n",
    "https://pandas.pydata.org/pandas-docs/stable/user_guide/groupby.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_aea_rgi.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_aea_rgi.groupby('RGIId')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hmmm.  Nothing happened.  Ah, we need a function to perform the aggregation over the grouped data!  How about taking the mean?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_aea_rgi.groupby('RGIId').mean().head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define a more sophisticated aggregation function\n",
    "\n",
    "A dictionary of fields and functions can be used to compute a set of summary statistics for relevant columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "agg_func = {'Name':'first',\n",
    "            'Area':'first', \n",
    "            'glas_z':['mean', 'median', 'std', nmad],\n",
    "            'glas_srtm_dhdt':['count','mean', 'median', 'std', nmad]}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_aea_rgi_agg = glas_gdf_aea_rgi.groupby('RGIId').agg(agg_func)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_aea_rgi_agg.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We don't need the multi-index here\n",
    "glas_gdf_aea_rgi_agg.columns = ['_'.join(col).rstrip('_') for col in glas_gdf_aea_rgi_agg.columns.values]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_aea_rgi_agg.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_aea_rgi_agg['glas_pt_density'] = glas_gdf_aea_rgi_agg['glas_srtm_dhdt_count']/glas_gdf_aea_rgi_agg['Area_first']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_aea_rgi_agg.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Wait a minute, what happened to our RGI polygon geometry?\n",
    "\n",
    "This was a casualty of our initial spatial join, as we preserved the Point geometry for each GLAS record, not the RGI geometry.\n",
    "\n",
    "Let's create a new GeoDataFrame, adding the original RGI geometry to the aggregated statistics.\n",
    "\n",
    "Since both DataFrames have the same Index (RGIId), Pandas will automatically join with corresponding records."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgi_gdf_conus_aea['geometry'].head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A bit of a hack to join on RGIId \n",
    "glas_gdf_aea_rgi_agg_gdf = gpd.GeoDataFrame(glas_gdf_aea_rgi_agg, geometry=rgi_gdf_conus_aea[['RGIId','geometry']].set_index('RGIId')['geometry'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_aea_rgi_agg_gdf.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Let's make some final maps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import some other useful mapping packages\n",
    "\n",
    "`matplotlib-scalebar` adds a dynamic scalebar to matplotlib axes\n",
    "\n",
    "`contextily` downloads and statically renders basemap tiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib_scalebar.scalebar import ScaleBar\n",
    "import contextily as ctx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_basemap(ax, crs, zoom='auto'):\n",
    "    ctx.add_basemap(ax=ax, crs=crs, url=ctx.sources.ST_TERRAIN, zoom=zoom)\n",
    "    #Create a scalebar object, with scaling factor of 1.0 px, since we're using projected coordinate system with unit 1 m\n",
    "    scalebar = ScaleBar(1.0)\n",
    "    #Add scalebar to axes\n",
    "    ax.add_artist(scalebar)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots()\n",
    "glas_gdf_aea.plot(ax=ax, column='glas_z', cmap='inferno', markersize=2, legend=True)\n",
    "#Add basemap, specifying crs keyword\n",
    "add_basemap(ax, crs=glas_gdf_aea.crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clim = (-5.0, 5.0)\n",
    "f, axa = plt.subplots(1,2, figsize=(10,4), sharex=True, sharey=True)\n",
    "\n",
    "rgi_gdf_conus_aea.plot(ax=axa[0], edgecolor='k', facecolor='w', lw=0.5, alpha=0.2);\n",
    "glas_gdf_aea_rgi_agg_gdf.plot(ax=axa[0],column='glas_srtm_dhdt_count', cmap='inferno', edgecolor='k', lw=0.5, legend=True)\n",
    "axa[0].set_title('Number of GLAS shots per RGI Polygon')\n",
    "add_basemap(axa[0], glas_gdf_aea.crs)\n",
    "\n",
    "glas_gdf_aea_rgi_agg_gdf.plot(ax=axa[1],column='glas_srtm_dhdt_median', cmap='RdBu', edgecolor='k', lw=0.5, vmin=clim[0], vmax=clim[1], legend=True)\n",
    "axa[1].set_title('Median SRTM to GLAS dh/dt (m/yr)')\n",
    "glas_gdf_aea_rgi.plot(ax=axa[1],column='glas_srtm_dhdt', cmap='RdBu', markersize=5, edgecolor='0.5', lw=0.5, vmin=clim[0], vmax=clim[1])\n",
    "states_gdf_aea.plot(ax=axa[1], edgecolor='k', facecolor='none', lw=0.5)\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Zoom to North Cascades"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Can we extract any additional insight from the reduced data?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Point count vs. Glacier Area\n",
    "\n",
    "*Do we see more valid GLAS points over bigger glaciers?*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots()\n",
    "ax.scatter(glas_gdf_aea_rgi_agg_gdf['Area_first'], glas_gdf_aea_rgi_agg_gdf['glas_srtm_dhdt_count'])\n",
    "ax.set_ylabel('Number of GLAS shots')\n",
    "ax.set_xlabel('Glacier Area (km^2)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Median dh/dt vs. Glacier Area\n",
    "\n",
    "Do we see less elevation change over bigger glaciers?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots()\n",
    "ax.scatter(glas_gdf_aea_rgi_agg_gdf['Area_first'], glas_gdf_aea_rgi_agg_gdf['glas_srtm_dhdt_median'])\n",
    "ax.set_xlabel('Glacier Area (km^2)')\n",
    "ax.set_ylabel('Median SRTM - GLAS dh/dt (m/yr)')\n",
    "ax.axhline(0, color='k');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Median dh/dt vs. Glacier Elevation\n",
    "\n",
    "*Do we see greater elevation loss at lower elevations?*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots()\n",
    "ax.scatter(glas_gdf_aea_rgi_agg_gdf['glas_z_median'], glas_gdf_aea_rgi_agg_gdf['glas_srtm_dhdt_median'])\n",
    "ax.set_xlabel('Median Glacier Elevation (m HAE)')\n",
    "ax.set_ylabel('Median SRTM - GLAS dh/dt (m/yr)')\n",
    "ax.axhline(0, color='k');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusion:\n",
    "\n",
    "N/A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save the final polygons to a GIS-ready file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The workflows in these Notebooks are intended to be fully reproducible, starting with raw data and producing all final output.  But sometimes you want to write out geospatial data for analysis in a GUI-based GIS (QGIS, ArcMap), or to share with colleagues who will use these tools."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check available output formats for geopandas\n",
    "* Use fiona to get a list of available file type drivers for output\n",
    "* Note: the 'r' means fiona/geopandas can read this file type, 'w' means it can write this file type, 'a' means it can append to an existing file.\n",
    "    * https://fiona.readthedocs.io/en/latest/manual.html#writing-vector-data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import fiona\n",
    "fiona.supported_drivers"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How to choose a format?\n",
    "* If you've taken a GIS class (or not), you've probably used shapefiles in the past.  Please stop.  The ESRI shapefile is a legacy format, though it is still widely used.\n",
    "* http://switchfromshapefile.org/\n",
    "* Better options these days are Geopackage (GPKG) when spatial index is required, and simple GeoJSON for other cases\n",
    "    * Both should be supported by your GIS (including QGIS, ArcGIS, etc)\n",
    "* Let's use geopackage for this exercise\n",
    "* Can use the Geopandas `to_file()` method to create this file\n",
    "    * Make sure you properly specify filename with extension and the `driver` option\n",
    "    * *Note: Writing out may take a minute, and may produce an intermediate '.gpkg-journal' file*\n",
    "        * Can see this in the file browser or terminal!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_fn='./conus_glas_gdf_aea_rgi_agg.gpkg'\n",
    "glas_gdf_aea_rgi_agg_gdf.to_file(out_fn, driver='GPKG')\n",
    "\n",
    "#out_fn='./conus_glas_gdf_aea_rgi_agg.geojson'\n",
    "#glas_gdf_aea_rgi_agg_gdf.to_file(out_fn, driver='GeoJSON')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ls -lh $out_fn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 🎉\n",
    "\n",
    "You can now directly load this gpkg file in any GIS, without defining a coordinate system. You can also load this file directly into geopandas in the future using the `read_file()` method, without having to repeat the processing above.\n",
    "\n",
    "### See for yourself!\n",
    "Try it! Right-click on file in the file browser to the left of the JupyterLab interface, then select Download and pick a location on your local computer (e.g., your Downloads folder). \n",
    "\n",
    "Then open this file in QGIS or ArcGIS on your local machine!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# So, what if we actually had a decent GLAS sample of CONUS glaciers?\n",
    "\n",
    "## Estimate glacier mass balance\n",
    "\n",
    "We could estimate mass balance for each glacier polygon using the mean dh/dt, glacier area, and a standard bulk density (e.g., 850 kg/m3 [Huss, 2013]).  Could then use geopandas to perform a secondary aggregation to compile statistics for polygons representing mountain ranges, river basins, etc.\n",
    "\n",
    "With a sparse sample, probably best to derive dh/dt vs. elevation curves, then combine with observed glacier hypsometry to estimate mass balance for glacier polygons.  This can work, but need to be careful about spatial variability."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Final thoughts\n",
    "\n",
    "GLAS was not the right tool for small, low-latitude CONUS glaciers.  We kind of knew this would be the case before we started, but hey, it was worth a look, and we learned some basic Python geospatial analysis skills.\n",
    "\n",
    "The concepts and approaches presented here can be applied to larger glaciers or ice caps, especially at higher latitudes.  One could modify to repeat for all RGI regions.\n",
    "\n",
    "We started with an existing csv of culled points here.  One could repeat with a similarly processed subset of ATL06 points using the workflows presented earlier this week.  This will provide a longer time period to evaluate noisy elevation measurements.  Replacing void-filled SRTM with another reference DEM is also needed (e.g., the timestamped USGS NED).\n",
    "\n",
    "Note that the core tools presented here have Dask integration (https://dask.org/) to allow you to chunk and process larger-than-memory datasets with minimal modification to the code."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other visualization packages for large point datasets\n",
    "* Datashader\n",
    "* ipyvolume\n",
    "* hvplot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Appendix A: Sampling a raster at points\n",
    "\n",
    "This is something that is surprisingly common, but may not be simple to implement.  Let's discuss a few options:\n",
    "1. Simple `rasterio` sampling with integer indices using nearest neighbor\n",
    "2. Statistics extracted for a circular window around each point location\n",
    "3. NumPy/SciPy interpolation routines"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this example, we will use a sample 90-m SRTM-GL1 DEM over WA state, but you could repeat with any DEM (e.g., ArcticDEM)\n",
    "\n",
    "Also, here is an `icepyx` notebook that includes an example of raster sampling: https://github.com/icesat2py/icepyx/blob/master/doc/examples/ICESat-2_DEM_comparison_Colombia_working.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Rasterio sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import requests\n",
    "import rasterio as rio\n",
    "from rasterio import plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fetch SRTM data\n",
    "* Let's use the excellent OpenTopography GlobalDEM REST API for this!\n",
    "* https://portal.opentopography.org/apidocs/#/operations/Public"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wa_bounds = (-124.733174, 45.543541, -116.915989, 49.002494)\n",
    "tile_bounds = (-125, 45, -120, 50)\n",
    "bounds = wa_bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "demtype = \"SRTMGL3\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "url=\"https://portal.opentopography.org/API/globaldem?demtype={}&west={}&south={}&east={}&north={}&outputFormat=GTiff\".format(demtype,*wa_bounds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "response = requests.get(url)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Open directly\n",
    "#src = rio.open(response.content)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Save to disk\n",
    "out_fn = 'wa_{}.tif'.format(demtype)\n",
    "if not os.path.exists(out_fn):\n",
    "    open(out_fn, 'wb').write(response.content)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Reproject to UTM\n",
    "proj_fn = os.path.splitext(out_fn)[0]+'_utm.tif'\n",
    "if not os.path.exists(proj_fn):\n",
    "    !gdalwarp -r cubic -t_srs EPSG:32610 -co COMPRESS=LZW -co TILED=YES -co BIGTIFF=IF_SAFER $out_fn $proj_fn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Open the file with rasterio\n",
    "\n",
    "https://rasterio.readthedocs.io/en/stable/"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "srtm_fn = proj_fn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#This is rasterio Dataset object\n",
    "srtm_src = rio.open(srtm_fn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "srtm_src.profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "srtm_src.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "srtm_extent = rio.plot.plotting_extent(srtm_src)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read as a NumPy Masked array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "srtm = srtm_src.read(1, masked=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "srtm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create a quick plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots()\n",
    "ax.imshow(srtm, extent=srtm_extent);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate shaded relief map\n",
    "* Many ways to do this, but we'll just use the `gdaldem` command line utility for simplicity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hs_fn = os.path.splitext(srtm_fn)[0]+'_hs.tif'\n",
    "if not os.path.exists(hs_fn):\n",
    "    !gdaldem hillshade $srtm_fn $hs_fn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "srtm_hs_src = rio.open(hs_fn)\n",
    "hs = srtm_hs_src.read(1, masked=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot color shaded relief map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots()\n",
    "ax.imshow(hs, cmap='gray', extent=rio.plot.plotting_extent(srtm_hs_src))\n",
    "ax.imshow(srtm, extent=srtm_extent, alpha=0.5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reproject GLAS points to match raster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_srtm = glas_gdf_aea.to_crs(srtm_src.crs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prepare the coordinate arrays to pass to rio `sample` function\n",
    "* The `sample` function expects a list of (x,y) tuples: https://rasterio.readthedocs.io/en/latest/api/rasterio.sample.html\n",
    "    * Need to create this from the `geometry` objects in your GeoDataFrame\n",
    "    * You want a list of the form [(x1,y1),(x2,y2),…]\n",
    "* Pass to `sample`\n",
    "* Note that the `sample` function returns a `generator` object, and it doesn't actually evaluate the call!\n",
    "* Can wrap this in a `np.array(list())` to evaluate, or use `np.fromiter()`\n",
    "* This operation may take ~10-20 seconds to complete"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_coord = [(pt.x, pt.y) for pt in glas_gdf_srtm.geometry]\n",
    "#glas_coord = np.vstack((glas_gdf_srtm.geometry.x.values, glas_gdf_srtm.geometry.y.values)).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sample with rasterio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_srtm_sample = srtm_src.sample(glas_coord)\n",
    "glas_srtm_sample"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### This is a generator, so we actually need to evaluate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_srtm_elev = np.fromiter(glas_srtm_sample, dtype=srtm.dtype)\n",
    "glas_srtm_elev"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Deal with nodata\n",
    "* Some of our GLAS points are located over areas where we don't have valid DEM pixels\n",
    "* These will be assigned the raster nodata value (-32768 in this case)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_srtm_elev_ma = np.ma.masked_equal(glas_srtm_elev, srtm_src.nodata)\n",
    "glas_srtm_elev_ma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add new column to the GeoDataFrame\n",
    "* Set masked values to `np.nan` (which requires a conversion to float)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_srtm['srtm_90m_z_rio'] = glas_srtm_elev_ma.astype(float).filled(np.nan)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_srtm.dropna().head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots()\n",
    "ax.imshow(hs, cmap='gray', extent=rio.plot.plotting_extent(srtm_hs_src))\n",
    "#ax.imshow(srtm, extent=srtm_extent, alpha=0.5);\n",
    "glas_gdf_srtm.dropna().plot('srtm_90m_z_rio', ax=ax, markersize=1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Note: the SRTM elevation values are height above the EGM96 geoid*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Notes on sampling coarse rasters or noisy rasters at integer pixel locations\n",
    "* The rasterio approach is efficient, but it uses a nearest neighbor algorithm to extract the elevation value for the grid cell that contains the point, regardless of where the point falls within the grid cell (center vs. corner)\n",
    "* But the DEM grid cells can be big (~90x90 m for the SRTM3 data), so if point is near the corner of a pixel on steep slope, the pixel value might not be representative.\n",
    "* A better approach is to use bilinear or bicubic sampling, to interpolate the elevation value at the point coordinate using pixel values within some neighborhood around the point, (e.g. 2x2 window for bilinear, 4x4 window for cubic)\n",
    "* Other approaches involve computing zonal stats within some radius of the point location (e.g., median elevation of pixels within 300 m of the point, using buffer to create polygons)\n",
    "    * https://www.earthdatascience.org/courses/earth-analytics-python/lidar-remote-sensing-uncertainty/extract-data-from-raster/\n",
    "    * https://pysal.org/scipy2019-intermediate-gds/deterministic/gds2-rasters.html#getting-values-at-cells\n",
    "    * https://github.com/dshean/pygeotools/blob/master/pygeotools/lib/geolib.py#L1019"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Local window sample\n",
    "\n",
    "https://github.com/dshean/demcoreg/blob/master/demcoreg/sample_raster_at_pts.py\n",
    "\n",
    "https://github.com/dshean/pygeotools/blob/master/pygeotools/lib/geolib.py#L1019"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Scipy ndimage: n-order polynomial\n",
    "* Good option for regular grids (i.e., DEMs)\n",
    "* Propagates nan, issues when DEM has missing data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.ndimage\n",
    "#Should dropna here\n",
    "yx = np.array([glas_gdf_srtm.geometry.y, glas_gdf_srtm.geometry.x])\n",
    "#Convert map coordinates to array coordinates (want float, not whole integer indices)\n",
    "#See functions here\n",
    "#https://rasterio.readthedocs.io/en/latest/api/rasterio.transform.html\n",
    "#Need to revisit\n",
    "#Use bilinear resampling here (order=1)\n",
    "ndimage_samp = scipy.ndimage.map_coordinates(srtm, yx, order=1, mode='nearest')\n",
    "ndimage_samp = np.ma.masked_equal(ndimage_samp, srtm_src.nodata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "srtm_src.transform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_srtm['srtm_90m_z_ndimage'] = ndimage_samp.astype(float).filled(np.nan)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "glas_gdf_srtm.dropna().head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Appendix B: Distortion in Polar Stereographic Projections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\"The [polar stereographic] projection is defined on the entire sphere, except at one point: the projection point. Where it is defined, the mapping is smooth and bijective. It is conformal, meaning that it preserves angles at which curves meet. It is neither isometric nor area-preserving: that is, it preserves neither distances nor the areas of figures.\" [Wikipedia]\n",
    "\n",
    "Many of you are probably using polar stereographic projections, likely the standard EPSG:3031 or EPSG:3413 projections for Antarctica and Greenland, respectively.  These are designed to minimize distortion near the latitude of true scale (71°S for Antarctica and 70°N for Greenland).  This means that area and distance distortion will increase as you move away from this latitude.  So areas or distances measured near the pole will not be equivalent to areas or distances measured on the Antarctic Peninsula.  The difference isn't huge, but is nontrivial (~1-2% for Antarctica I believe).  See https://nsidc.org/data/polar-stereo/ps_grids.html for more details.\n",
    "\n",
    "The figures here might help illustrate: https://en.wikipedia.org/wiki/Stereographic_projection#Properties\n",
    "\n",
    "Let's try to illustrate this nuanced, but often overlooked issue."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tissot's Indicatrix for polar stereographic projection\n",
    "\n",
    "Let's use the classic Tissot's Indicatrix plots to show map distortion.\n",
    "\n",
    "https://en.wikipedia.org/wiki/Tissot%27s_indicatrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cartopy\n",
    "import cartopy.crs as ccrs\n",
    "\n",
    "#Cartopy implementation of EPSG:3031 and EPSG:3413\n",
    "crs_3031 = ccrs.SouthPolarStereo(true_scale_latitude=-71)\n",
    "crs_3413 = ccrs.NorthPolarStereo(central_longitude=-45, true_scale_latitude=70)\n",
    "\n",
    "#Circle locations\n",
    "lons = range(-180, 180, 30)\n",
    "lats = range(-90, 91, 10)\n",
    "#Radius of circles\n",
    "rad = 400"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10,5))\n",
    "ax1 = plt.subplot(1, 2, 1, projection=crs_3413)\n",
    "ax1.coastlines()\n",
    "ax1.gridlines(ylocs=[70,],color='r')\n",
    "ax1.tissot(facecolor='orange', edgecolor='0.5', alpha=0.4, rad_km=rad, lons=lons, lats=lats)\n",
    "ax1.set_title('EPSG:3413')\n",
    "\n",
    "ax2 = plt.subplot(1, 2, 2, projection=crs_3031)\n",
    "ax2.coastlines()\n",
    "ax2.gridlines(ylocs=[-71,],color='r')\n",
    "ax2.tissot(facecolor='orange', edgecolor='0.5', alpha=0.4, rad_km=rad, lons=lons, lats=lats)\n",
    "ax2.set_title('EPSG:3031')\n",
    "ax2.set_extent([-180, 180, -90, -50], ccrs.PlateCarree())\n",
    "\n",
    "#There is a bug in cartopy extent when central_longitude is not 0\n",
    "#Get extent in projected crs (x0, x1, y0, y1)\n",
    "extent_3413 = ax2.get_extent()\n",
    "ax1.set_extent(extent_3413, crs_3413)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OK, cool plot bro.  But why does this matter?\n",
    "\n",
    "Note the size of the circles in the corners and over the pole, relative to the circles near the latitude of true scale (red line).  While it is unlikely that you'll use this projection to look at mid-latitude regions, you can see the difference in area distortion over the ~20° of latitude between North and South Greenland.\n",
    "\n",
    "Say you generated amazing elevation difference grids for all of the Greenland ice sheet using ICESat and ICESat-2 crossovers, and you used the EPSG:3413 projection.  Say the grids have a posting (grid cell size) of 120 m.  Imagine didividing the above plots into grid cells.  You'll end up with more grid cells over the circles at lower latitudes (e.g., South Greenland) and fewer grid cells over circles at high latitudes (e.g., North Greenland).  \n",
    "\n",
    "Let's put it another way.  You could try to compute volume change of the Greenland ice sheet by summing dh/dt values for all grid cells in each catchment.  Let's assume all grid cells have the same -1.0 m/yr value.  The integrated estimates for the catchments in South Greenland will have more grid cells due to the projection, resulting in larger apparent negative volume change!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Let's compare with an equal-area projection\n",
    "Lambert Azimuthal Equal-Area is not a bad choice for the polar regions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Cartopy implementation of EPSG:3031 and EPSG:3413\n",
    "#crs_slaea = ccrs.LambertAzimuthalEqualArea(central_longitude=0.0, central_latitude=-90.0)\n",
    "crs_nlaea = ccrs.LambertAzimuthalEqualArea(central_longitude=-45.0, central_latitude=90.0)\n",
    "\n",
    "#Specify locations of Tissot's Indicatrix\n",
    "lons = range(-180, 180, 30)\n",
    "lats = range(-90, 91, 10)\n",
    "rad = 400"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Add a 100 km grid\n",
    "dx = 100000\n",
    "xgrid = np.arange(extent_3413[0], extent_3413[1]+dx, dx)\n",
    "ygrid = np.arange(extent_3413[2], extent_3413[3]+dx, dx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10,5))\n",
    "ax1 = plt.subplot(1, 2, 1, projection=crs_3413)\n",
    "ax1.coastlines()\n",
    "ax1.tissot(facecolor='orange', edgecolor='0.5', alpha=0.4, rad_km=rad, lons=lons, lats=lats)\n",
    "ax1.set_title('EPSG:3413')\n",
    "#ax1.gridlines(crs=ccrs.PlateCarree(), ylocs=lats)\n",
    "ax1.gridlines(crs=crs_3413, xlocs=xgrid, ylocs=ygrid)\n",
    "ax1.gridlines(ylocs=[70,],color='r')\n",
    "\n",
    "ax2 = plt.subplot(1, 2, 2, projection=crs_nlaea)\n",
    "ax2.coastlines()\n",
    "ax2.tissot(facecolor='orange', edgecolor='0.5', alpha=0.4, rad_km=rad, lons=lons, lats=lats)\n",
    "ax2.set_title('Lambert Azimuthal Equal-Area')\n",
    "ax2.set_extent([-180, 180, 50, 90], ccrs.PlateCarree())\n",
    "#ax2.gridlines(crs=ccrs.PlateCarree(), ylocs=lats)\n",
    "ax2.gridlines(crs=crs_nlaea, xlocs=xgrid, ylocs=ygrid)\n",
    "\n",
    "#ax1.set_extent(extent_3413, crs_3413)\n",
    "#ax2.set_extent(extent_3413, crs_3413)\n",
    "\n",
    "gr_extent = [-1.4E6,1.4E6,-3.7E6,-5.3E5]\n",
    "ax1.set_extent(gr_extent, crs_3413)\n",
    "ax2.set_extent(gr_extent, crs_3413)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not the distorted shape of the circles near Southern Greenland, though they should have identical area."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bonus: the \"Flat Earth\" Projection: Azimuthal Equidistant"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Azimuthal Equidistant](https://upload.wikimedia.org/wikipedia/commons/2/2f/Flat_earth.png)\n",
    "\n",
    ">An azimuthal equidistant projection of the entire spherical Earth. A rendered picture of the Flat Earth model. The white around the outside of the globe is thought to be an 'Ice Wall', preventing people from falling off the surface of the earth. [Wikipedia]\n",
    "\n",
    "So, those of you who have been to the South Pole, can you help me understand how this works?"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}