{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Sliderule data I/O\n",
    "\n",
    "This notebook will highlight analyzing various coincident elevation measurements. We will find regions with and use sliderule to retrieve ICESat-2 and GEDI point elevation measurements.\n",
    "\n",
    "```{note}\n",
    "Keep in mind, these measurements are from different sensor types, close in time, but not at exactly the same time. All measurements also have uncertainties, so we do not expect perfect agreement.\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import coincident\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "# For testing\n",
    "# import sliderule\n",
    "# sliderule.init(url='slideruleearth.io', verbose=True)\n",
    "\n",
    "from mpl_toolkits.axes_grid1.inset_locator import inset_axes\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",
    "gf_lidar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Search secondary datasets\n",
    "\n",
    "Provide a list that will be searched in order. The list contains tuples of dataset aliases and the temporal pad in days to search before the primary dataset start and end dates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "secondary_datasets = [\n",
    "    (\"gedi\", 40),  # +/- 40 days from lidar\n",
    "    (\"icesat-2\", 60),  # +/- 60 days from lidar\n",
    "]\n",
    "\n",
    "gf_gedi, gf_is2 = coincident.search.cascading_search(\n",
    "    gf_lidar,\n",
    "    secondary_datasets,\n",
    "    min_overlap_area=30,  # km^2\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get ICESat-2 ATL06 Data\n",
    "\n",
    "We've identified 7 granules of icesat-2 data to examine, but there is no need to work with the entire granule, which spans a huge geographic extent. Instead we'll retrieve a subset of elevation values in the area of interest for each granule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hone in on single granule as a pandas series\n",
    "i = 0\n",
    "granule_gdf = gf_is2.iloc[[i]]\n",
    "granule = gf_is2.iloc[i]\n",
    "granule_gdf.plot()  # needs to be a geodataframe\n",
    "plt.title(granule.id);\n",
    "# The icesat-2 track trends N-S, and is crossed by multiple GEDI tracks, resulting in the crosshatched appearance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_is2 = coincident.io.sliderule.subset_atl06(\n",
    "    granule_gdf,\n",
    "    include_worldcover=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot this data and overlay footprint\n",
    "fig, ax = plt.subplots(figsize=(8, 10))\n",
    "plt.scatter(\n",
    "    x=data_is2.geometry.x,\n",
    "    y=data_is2.geometry.y,\n",
    "    c=data_is2.h_li,\n",
    "    s=10,\n",
    ")\n",
    "granule_gdf.dissolve().boundary.plot(ax=ax, color=\"magenta\")\n",
    "cb = plt.colorbar()\n",
    "cb.set_label(\"elevation_hr (m)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```{note}\n",
    "`sliderule` gets data from the envelope of multipolygon, not only data within intersecting GEDI tracks. Along-track gaps are where there is missing data.\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here, we specify aoi=granule_gdf to only return GEDI data\n",
    "# in our icesat-2 track of interest\n",
    "data_gedi = coincident.io.sliderule.subset_gedi02a(\n",
    "    gf_gedi,\n",
    "    aoi=granule_gdf,\n",
    "    include_worldcover=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data_gedi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot this data and overlay footprint\n",
    "fig, ax = plt.subplots(figsize=(8, 10))\n",
    "plt.scatter(\n",
    "    x=data_gedi.geometry.x,\n",
    "    y=data_gedi.geometry.y,\n",
    "    c=data_gedi.elevation_hr,\n",
    "    s=10,\n",
    ")\n",
    "granule_gdf.dissolve().boundary.plot(ax=ax, color=\"magenta\")\n",
    "cb = plt.colorbar()\n",
    "cb.set_label(\"elevation_hr (m)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**NOTE** Like ICESat-2, it's possible to have points falling outside the estimated footprints"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Spatial join: nearest points\n",
    "\n",
    "NOTE: we will not worry about the difference in time of acquisition between adjacent points for now"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "utm_crs = granule_gdf.estimate_utm_crs()\n",
    "data_is2_utm = data_is2.to_crs(utm_crs)\n",
    "data_gedi_utm = data_gedi.to_crs(utm_crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# find nearest IS2 point to each GEDI\n",
    "close_gedi = data_gedi_utm.sjoin_nearest(\n",
    "    data_is2_utm,\n",
    "    how=\"left\",\n",
    "    max_distance=100,  # at most 100m apart\n",
    "    distance_col=\"distances\",\n",
    ")\n",
    "# index is GEDI subset (with _left added), with _right appended to is2 columns + Distances in meters\n",
    "close_gedi = close_gedi[close_gedi[\"distances\"].notna()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "plt.scatter(close_gedi.h_li, close_gedi.elevation_lm)\n",
    "ax.axline((0, 0), slope=1, color=\"k\", transform=ax.transAxes)\n",
    "plt.xlabel(\"ATL06 elevation (m)\")\n",
    "plt.ylabel(\"GEDI L2A elevation_lm (m)\")\n",
    "plt.title(f\"{len(close_gedi)} points within 100m\")\n",
    "\n",
    "\n",
    "# Add statistics and inset histogram to plot\n",
    "inset_ax = inset_axes(ax, width=\"40%\", height=\"30%\", borderpad=2, loc=\"lower right\")\n",
    "\n",
    "\n",
    "residuals = close_gedi.h_li - close_gedi.elevation_lm\n",
    "mean_residual = np.mean(residuals)\n",
    "median_residual = np.median(residuals)\n",
    "std_residual = np.std(residuals)\n",
    "textstr = \"\\n\".join(\n",
    "    (\n",
    "        f\"Mean: {mean_residual:.2f}\",\n",
    "        f\"Median: {median_residual:.2f}\",\n",
    "        f\"Std: {std_residual:.2f}\",\n",
    "    )\n",
    ")\n",
    "props = dict(boxstyle=\"round\", facecolor=\"wheat\", alpha=0.5)\n",
    "ax.text(\n",
    "    0.05,\n",
    "    0.95,\n",
    "    textstr,\n",
    "    transform=ax.transAxes,\n",
    "    fontsize=10,\n",
    "    verticalalignment=\"top\",\n",
    "    bbox=props,\n",
    ")\n",
    "\n",
    "_ = inset_ax.hist(residuals, bins=50, range=(-20, 20), color=\"gray\")\n",
    "inset_ax.axvline(0, color=\"k\", linestyle=\"-\", linewidth=1);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sample 3DEP\n",
    "\n",
    "Sample 3DEP 1m DEM at the subset of GEDI Points\n",
    "\n",
    "```{note}\n",
    "SlideRule returns all elevation data in EPSG:7912\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "close_gedi_3dep = coincident.io.sliderule.sample_3dep(\n",
    "    close_gedi,\n",
    "    # Restrict to only DEMs derived from specific WESM LiDAR project\n",
    "    project_name=gf_lidar[\"project\"].iloc[0],\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add 3DEP elevation values to existing dataframe\n",
    "close_gedi[\"elevation_3dep\"] = close_gedi_3dep.value.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Geometries of nearest IS2 points\n",
    "close_is2 = data_is2_utm.loc[close_gedi.time_right]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sample 3DEP at the subset of GEDI Points\n",
    "close_is2_3dep = coincident.io.sliderule.sample_3dep(\n",
    "    close_is2,\n",
    "    project_name=gf_lidar[\"project\"].iloc[0],\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "close_is2[\"elevation_3dep\"] = close_is2_3dep.value.values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.axvline(0, color=\"k\", linestyle=\"dashed\", linewidth=1)\n",
    "\n",
    "diff_3dep = close_is2.h_li.values - close_is2.elevation_3dep\n",
    "label = f\"ICESat-2 (median={np.nanmedian(diff_3dep):.2f}, mean={np.nanmean(diff_3dep):.2f}, std={np.nanstd(diff_3dep):.2f}, n={len(diff_3dep)})\"\n",
    "plt.hist(diff_3dep, range=[-5, 5], bins=100, color=\"m\", alpha=0.5, label=label)\n",
    "\n",
    "diff_3dep = close_gedi.elevation_lm.values - close_gedi.elevation_3dep\n",
    "label = f\"GEDI2A (median={np.nanmedian(diff_3dep):.2f}, mean={np.nanmean(diff_3dep):.2f}, std={np.nanstd(diff_3dep):.2f}, n={len(diff_3dep)})\"\n",
    "plt.hist(diff_3dep, range=[-5, 5], bins=100, color=\"c\", alpha=0.5, label=label)\n",
    "\n",
    "plt.xlabel(\"elevation difference (m)\")\n",
    "plt.title(\"Altimeter Elevation - 3DEP 1m DEM\")\n",
    "plt.xlim(-5, 5)\n",
    "plt.legend(loc=\"upper left\");"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "docs",
   "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}