{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Contextual data\n",
    "\n",
    "Once you've identified areas of interest where multiple datasets intersect, you can pull additional data to provide further context. For example:\n",
    "\n",
    "1. landcover \n",
    "2. global elevation data "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import coincident\n",
    "import matplotlib.pyplot as plt\n",
    "import xarray as xr\n",
    "import numpy as np\n",
    "import rasterio\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Identify a primary dataset\n",
    "\n",
    "Start by loading a full resolution polygon of a 3DEP LiDAR workunit which has a known start_datetime and end_datatime:"
   ]
  },
  {
   "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": [
    "search_aoi = gf_lidar.simplify(0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Search for coincident contextual data\n",
    "\n",
    "Coincident provides two convenience functions to load datasets resulting from a given search. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gf_wc = coincident.search.search(\n",
    "    dataset=\"worldcover\",\n",
    "    intersects=search_aoi,\n",
    "    # worldcover is just 2020 an 2021, so pick one\n",
    "    datetime=[\"2020\"],\n",
    ")  # Asset of interest = 'map'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# STAC metadata always has a \"stac_version\" column\n",
    "gf_wc.iloc[0].stac_version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gf_cop30 = coincident.search.search(\n",
    "    dataset=\"cop30\",\n",
    "    intersects=search_aoi,\n",
    ")  # Asset of interest = 'data'\n",
    "gf_cop30.iloc[0].stac_version"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### STAC search results to datacube\n",
    "\n",
    "If the results have [STAC-formatted metadata](https://stacspec.org/en). We can take advantage of the excellent [odc.stac](https://odc-stac.readthedocs.io/) tool to load datacubes. Please refer to odc.stac documentation for all the configuration options (e.g. setting resolution our output CRS, etc.)\n",
    "\n",
    "By default this uses the [dask](https://www.dask.org/) parallel computing library to intelligently load only metadata initially and defer reading values until computations or visualizations are performed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Copernicus DEM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = coincident.io.xarray.to_dataset(\n",
    "    gf_cop30,\n",
    "    aoi=search_aoi,\n",
    "    # chunks=dict(x=2048, y=2048), # manual chunks\n",
    "    resolution=0.00081,  # ~90m\n",
    "    mask=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# By default, these are dask arrays\n",
    "ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# You might want to rename the data variable\n",
    "ds = ds.rename(data=\"elevation\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The total size of this dataset is only 14MB, so for faster computations, load it into memory\n",
    "print(ds.nbytes / 1e6)\n",
    "ds = ds.compute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.elevation.isel(time=0).plot.imshow();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### ESA Worldcover"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Same with LandCover\n",
    "dswc = coincident.io.xarray.to_dataset(\n",
    "    gf_wc,\n",
    "    bands=[\"map\"],\n",
    "    aoi=search_aoi,\n",
    "    mask=True,\n",
    "    # resolution=0.00027, #~30m\n",
    "    resolution=0.00081,  # ~90m\n",
    ")\n",
    "dswc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dswc = dswc.rename(map=\"landcover\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dswc = dswc.compute()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For landcover there is a convenicence function for a nice categorical colormap\n",
    "ax = coincident.plot.plot_esa_worldcover(dswc)\n",
    "ax.set_title(\"ESA WorldCover\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load gridded elevation in a consistent CRS\n",
    "\n",
    "`coincident` also has a convenience function for loading gridded elevation datasets in a consistent CRS. In order to facilitate comparison with modern altimetry datasets (ICESat-2, GEDI), we convert elevation data on-the-fly to [EPSG:7912](https://spatialreference.org/ref/epsg/7912/). 3D Coordinate Reference Frames are a complex topic, so please see this resource for more detail on how these conversions are done: https://uw-cryo.github.io/3D_CRS_Transformation_Resources/ \n",
    "\n",
    "```{note}\n",
    "Currently gridded DEMs are retrieved from [OpenTopography](https://opentopography.org/) hosted in AWS us-west-2.\n",
    "```\n",
    "\n",
    "```{warning}\n",
    " The output of `load_dem_7912()` uses the native resolution and grid of the input dataset, and data is immediately read into memory. So this method is better suited to small AOIs.\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Start with a small AOI:\n",
    "from shapely.geometry import box\n",
    "import geopandas as gpd\n",
    "\n",
    "aoi = gpd.GeoDataFrame(\n",
    "    geometry=[box(-106.812163, 38.40825, -106.396812, 39.049796)], crs=\"EPSG:4326\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "da_cop = coincident.io.xarray.load_dem_7912(\"cop30\", aoi=aoi)\n",
    "da_cop"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# NASA DEM and COP30 DEM are both 30m but will not necessarily have same coordinates!\n",
    "da_nasa = coincident.io.xarray.load_dem_7912(\"nasadem\", aoi=aoi)\n",
    "da_nasa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We check that coordinates are sufficiently close before force-alignment.\n",
    "# If sufficiently different, it would be better to use rioxarray.reproject_match!\n",
    "xr.testing.assert_allclose(da_cop.x, da_nasa.x)  # rtol=1e-05, atol=1e-08 defaults\n",
    "xr.testing.assert_allclose(da_cop.y, da_nasa.y)\n",
    "da_nasa = da_nasa.assign_coords(x=da_cop.x, y=da_cop.y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diff = da_cop - da_nasa\n",
    "median = diff.median()\n",
    "mean = diff.mean()\n",
    "std = diff.std()\n",
    "\n",
    "fig = plt.figure(layout=\"constrained\", figsize=(4, 8))  # figsize=(8.5, 11)\n",
    "ax = fig.add_gridspec(top=0.75).subplots()\n",
    "ax_histx = ax.inset_axes([0, -0.35, 1, 0.25])\n",
    "axes = [ax, ax_histx]\n",
    "diff.plot.imshow(\n",
    "    ax=axes[0], robust=True, add_colorbar=False\n",
    ")  # , cbar_kwargs={'label': \"\"})\n",
    "n, bins, ax = diff.plot.hist(ax=axes[1], bins=100, range=(-20, 20), color=\"gray\")\n",
    "approx_mode = bins[np.argmax(n)]\n",
    "axes[0].set_title(\"COP30 - NASADEM\")\n",
    "axes[0].set_aspect(aspect=1 / np.cos(np.deg2rad(38.7)))\n",
    "axes[1].axvline(0, color=\"k\")\n",
    "axes[1].axvline(median, label=f\"median={median:.2f}\", color=\"cyan\", lw=1)\n",
    "axes[1].axvline(mean, label=f\"mode={mean:.2f}\", color=\"magenta\", lw=1)\n",
    "axes[1].axvline(approx_mode, label=f\"mode={approx_mode:.2f}\", color=\"yellow\", lw=1)\n",
    "axes[1].set_xlabel(\"Elevation difference (m)\")\n",
    "axes[1].legend()\n",
    "axes[1].set_title(\"\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{note}\n",
    "We see that COP30 elevation values are approximately 1m greater than NASADEM values for this area. Such a difference is within the stated accuracies of each gridded dataset, and also COP30 is derived from X-band TanDEM-X observations with nominal time of 2021-04-22 whereas NASADEM is primarily based on C-band SRTM collected 2000-02-20.\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load 3DEP 10m data\n",
    "da_3dep = coincident.io.xarray.load_dem_7912(\"3dep\", aoi=aoi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To compare to 3dep, we refine cop30 from 30m to 10m using bilinear resampling\n",
    "da_cop_r = da_cop.rio.reproject_match(\n",
    "    da_3dep, resampling=rasterio.enums.Resampling.bilinear\n",
    ")\n",
    "da_cop_r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "diff = da_3dep - da_cop_r\n",
    "median = diff.median()\n",
    "mean = diff.mean()\n",
    "std = diff.std()\n",
    "\n",
    "fig = plt.figure(layout=\"constrained\", figsize=(4, 8))  # figsize=(8.5, 11)\n",
    "ax = fig.add_gridspec(top=0.75).subplots()\n",
    "ax_histx = ax.inset_axes([0, -0.35, 1, 0.25])\n",
    "axes = [ax, ax_histx]\n",
    "diff.plot.imshow(\n",
    "    ax=axes[0], robust=True, add_colorbar=False\n",
    ")  # , cbar_kwargs={'label': \"\"})\n",
    "n, bins, ax = diff.plot.hist(ax=axes[1], bins=100, range=(-20, 20), color=\"gray\")\n",
    "approx_mode = bins[np.argmax(n)]\n",
    "axes[0].set_title(\"3DEP - COP30\")\n",
    "axes[0].set_aspect(aspect=1 / np.cos(np.deg2rad(38.7)))\n",
    "axes[1].axvline(0, color=\"k\")\n",
    "axes[1].axvline(median, label=f\"median={median:.2f}\", color=\"cyan\", lw=1)\n",
    "axes[1].axvline(mean, label=f\"mode={mean:.2f}\", color=\"magenta\", lw=1)\n",
    "axes[1].axvline(approx_mode, label=f\"mode={approx_mode:.2f}\", color=\"yellow\", lw=1)\n",
    "axes[1].set_xlabel(\"Elevation difference (m)\")\n",
    "axes[1].legend()\n",
    "axes[1].set_title(\"\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{note}\n",
    "There is a clear spatial and terrain dependence for the residuals of 10m 3DEP LiDAR compared to COP30 elevation values for this area! \n",
    "```"
   ]
  }
 ],
 "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
}