{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparing Elevation Measurements\n",
    "\n",
    "Once you have an area of interest with multiple overlapping DEMs and/or laser altimetry measurements, you can quickly compare elevation values. `coincident` supports a main plotting panel for elevation comparisons, as well as dynamic boxplots and histograms for elevation differences over various terrain characteristics and land cover types. \n",
    "\n",
    "This notebook highlights the following functions:\n",
    "\n",
    "* coincident.plot.compare_dems\n",
    "* coincident.plot.boxplot_slope\n",
    "* coincident.plot.boxplot_elevation\n",
    "* coincident.plot.boxplot_aspect\n",
    "* coincident.plot.hist_esa\n",
    "\n",
    "```{warning}\n",
    "Plotting routines are currently only tested for comparing datasets in EPSG:7912 (ITRF 2014)\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import coincident\n",
    "import geopandas as gpd\n",
    "from shapely.geometry import box\n",
    "import rioxarray as rxr\n",
    "import xarray as xr\n",
    "import xdem\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from rasterio.enums import Resampling\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load reference DEM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "workunit = \"CO_WestCentral_2019\"\n",
    "df_wesm = coincident.search.wesm.read_wesm_csv()\n",
    "gf_lidar = coincident.search.wesm.load_by_fid(\n",
    "    df_wesm[df_wesm.workunit == workunit].index\n",
    ")\n",
    "\n",
    "gf_lidar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We will examine the 'sourcedem' 1m\n",
    "# NOTE: it's important to note the *vertical* CRS of the data\n",
    "gf_lidar.iloc[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NOTE: since we're going to compare to satellite altimetry, let's convert first to EPSG:7912 !\n",
    "# GDAL / PROJ will use the following transform by default\n",
    "!projinfo -s EPSG:6350+5703 -t EPSG:7912 -q -o proj"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's grab the vendor DEM products for two arbitrary tiles of this flight to look at elevation differences between the aerial LIDAR derved DEM, Copernicus GLO-30 DEM, ICESat-2 elevations, and GEDI elevations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We will rely on GDAL to do this reprojection\n",
    "# gdalvsi = \"/vsicurl?empty_dir=yes&url=\"\n",
    "gdalvsi = \"/vsicurl/\"\n",
    "url_tile_1 = f\"{gdalvsi}https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/OPR/Projects/CO_WestCentral_2019_A19/CO_WestCentral_2019/TIFF/USGS_OPR_CO_WestCentral_2019_A19_be_w1032n1780.tif\"\n",
    "url_tile_2 = f\"{gdalvsi}https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/OPR/Projects/CO_WestCentral_2019_A19/CO_WestCentral_2019/TIFF/USGS_OPR_CO_WestCentral_2019_A19_be_w1032n1779.tif\"\n",
    "\n",
    "# First build a mosaic in native 3D CRS\n",
    "!gdalbuildvrt -q -a_srs EPSG:6350+5703 /tmp/mosaic.vrt '{url_tile_1}' '{url_tile_2}'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Then reproject\n",
    "!GDAL_DISABLE_READDIR_ON_OPEN=EMPTY_DIR gdalwarp -q -overwrite -t_srs EPSG:7912 -r bilinear /tmp/mosaic.vrt /tmp/mosaic_7912.tif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open our 'reference' DEM\n",
    "da_lidar = rxr.open_rasterio(\"/tmp/mosaic_7912.tif\", masked=True).squeeze()\n",
    "da_lidar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot a single DEM\n",
    "ax = coincident.plot.plot_dem(da_lidar, title=\"USGS_OPR_CO_WestCentral_2019_A19\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load secondary DEMs\n",
    "\n",
    "```{note}\n",
    "Currently coincident facilitates loading select global DEMs in EPSG:7912, if your dataset is not in that CRS be sure to first reproject your data\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# search geodataframe to grab the COP30 DEM\n",
    "bbox_tile = da_lidar.rio.bounds()\n",
    "poly_tile = box(*bbox_tile)\n",
    "gf_tile = gpd.GeoDataFrame(geometry=[poly_tile], crs=da_lidar.rio.crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NOTE: this is 30m data, but we are going to resample it to 1m to match our lidar\n",
    "da_cop = coincident.io.xarray.load_dem_7912(\"cop30\", aoi=gf_tile)\n",
    "da_cop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_cop_r = (\n",
    "    da_cop.rio.reproject_match(\n",
    "        da_lidar,\n",
    "        resampling=Resampling.bilinear,\n",
    "    )\n",
    "    .where(da_lidar.notnull())\n",
    "    .to_dataset(name=\"elevation\")\n",
    ")\n",
    "\n",
    "ds_cop_r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create hillshade variables for plot backgrounds\n",
    "# This function expects Datasets\n",
    "ds_lidar = da_lidar.to_dataset(name=\"elevation\")\n",
    "\n",
    "ds_lidar[\"hillshade\"] = coincident.plot.hillshade(ds_lidar.elevation)\n",
    "ds_cop_r[\"hillshade\"] = coincident.plot.hillshade(ds_cop_r.elevation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot a single DEM with hillshade\n",
    "ax = coincident.plot.plot_dem(\n",
    "    ds_cop_r[\"elevation\"],\n",
    "    title=\"Resampled Copernicus DEM\",\n",
    "    da_hillshade=ds_lidar[\"hillshade\"],\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get GEDI\n",
    "gf_gedi = coincident.search.search(\n",
    "    dataset=\"gedi\",\n",
    "    intersects=gf_tile,\n",
    ")\n",
    "data_gedi = coincident.io.sliderule.subset_gedi02a(\n",
    "    gf_gedi, gf_tile, include_worldcover=True\n",
    ")\n",
    "# Get ICESAT-2\n",
    "gf_is2 = coincident.search.search(\n",
    "    dataset=\"icesat-2\",\n",
    "    intersects=gf_tile,\n",
    ")\n",
    "data_is2 = coincident.io.sliderule.subset_atl06(\n",
    "    gf_is2, gf_tile, include_worldcover=True\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NOTE: this works because ICESat-2 and GEDI are EPSG:7912\n",
    "coincident.plot.plot_altimeter_points(\n",
    "    data_is2, \"h_li\", da_hillshade=ds_lidar[\"hillshade\"], title=\"ICESat-2 ATL06 points\"\n",
    ");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Main Elevation Comparison Panel\n",
    "\n",
    "'coincident' supports a main `compare_dems` function that allows you to compare multiple DEMs and altimetry points in a single panel. This function dynamically adjusts based on the number of datasets provided.\n",
    "\n",
    "```{note}\n",
    "`compare_dems` assumes that all DEMs are in the same coordinate reference system (CRS) and aligned. It also assumes that each DEM has an 'elevation' variable. If add_hillshade=True, it assumes the first DEM passed has a 'hillshade' variable, where the 'hillshade' variable can be calculated with `coincident.plot.create_hillshade` as seen below. A minimum of two DEMs must be provided and a maximum of 5 total datasets (datasets being DEMs and GeoDataFrames of altimetry points) can be provided.\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# list of DEMs to compare, where the first DEM in the list is the 'source' DEM\n",
    "dem_list = [ds_lidar, ds_cop_r]\n",
    "# dictionary of GeoDataFrames with altimetry points to compare\n",
    "# where the key is the name/plot title of the dataset and the value is a tuple of the respective\n",
    "# GeoDataFrame and the column name of the elevation values of interest\n",
    "dicts_gds = {\"ICESat-2\": (data_is2, \"h_li\"), \"GEDI\": (data_gedi, \"elevation_lm\")}\n",
    "\n",
    "axd = coincident.plot.compare_dems(\n",
    "    dem_list,\n",
    "    dicts_gds,\n",
    "    # ds_wc=ds_wc,\n",
    "    add_hillshade=True,\n",
    "    diff_clim=(-5, 5),  # difference colormap limits\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dynamic Elevation Difference Boxplots Based on Terrain Variables (slope, elevation, aspect)\n",
    "\n",
    "`coincident` allows users to dynamically plot elevation differences based on grouped data such as slope values. In these functions (that are wrappers for the generalized coincident.plot.boxplot_terrain_diff function), the user provides a reference DEM and one other DEM or GeoDataFrame with elevation values to compare to. These function assumes that the reference DEM has the respective terrain variable (slope, aspect, etc.) of interest.\n",
    "\n",
    "```{note}\n",
    "Boxplots will only be generated with groups of value counts greater than 30. Dynamic ylims will be fit to the IQRs of the boxplots.\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# THese routines assume rectilinear data, so we first convert to UTM\n",
    "crs_utm = ds_lidar.rio.estimate_utm_crs()\n",
    "ds_lidar_utm = ds_lidar.rio.reproject(crs_utm)\n",
    "ds_cop_utm = ds_cop_r.rio.reproject(crs_utm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope = xdem.terrain.slope(\n",
    "    ds_lidar_utm.elevation.values, resolution=ds_lidar_utm.rio.resolution()[0]\n",
    ")\n",
    "ds_lidar_utm[\"slope\"] = ((\"y\", \"x\"), np.squeeze(slope))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_, ax = plt.subplots()\n",
    "ds_lidar_utm.slope.plot.imshow(ax=ax)\n",
    "ax.set_title(\"USGS_OPR_CO_WestCentral_2019_A19\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax_boxplot_dems = coincident.plot.boxplot_slope([ds_lidar_utm, ds_cop_utm])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we also support elevation differences based on elevation bins\n",
    "# grouped on the reference DEM's elevation values\n",
    "ax_elev = coincident.plot.boxplot_elevation(\n",
    "    [ds_lidar_utm, ds_cop_utm], elevation_bins=np.arange(1700, 2000, 20)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# lastly for the boxplots, let's look at aspect\n",
    "aspect = xdem.terrain.aspect(ds_lidar_utm.elevation.values)\n",
    "ds_lidar_utm[\"aspect\"] = ((\"y\", \"x\"), np.squeeze(aspect))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "_, ax = plt.subplots()\n",
    "ds_lidar_utm.aspect.plot.imshow(ax=ax)\n",
    "ax.set_title(\"USGS_OPR_CO_WestCentral_2019_A19\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax_aspect = coincident.plot.boxplot_aspect([ds_lidar_utm, ds_cop_utm])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Histograms of Elevation Differences over Land Cover\n",
    "\n",
    "Plot elevation differences between DEMs or point data, grouped by ESA World Cover 2020 land cover class\n",
    "\n",
    "```{note}\n",
    "Histograms will only be generated with groups of value counts greater than a user-defined threshold (default is 30).\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax_lulc_cop30 = coincident.plot.hist_esa([ds_lidar, ds_cop_r])"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "dev",
   "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.12.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}