{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ddce2704-937f-46be-a0c8-fd402cc50037",
   "metadata": {},
   "source": [
    "## ICESat-2 ATL03 SlideRule Demo\n",
    "\n",
    "Plot ATL03 data with different classifications for a region over the Grand Mesa, CO region \n",
    "\n",
    "- [ATL08 Land and Vegetation Height product](https://nsidc.org/data/atl08) photon classification\n",
    "- Experimental YAPC (Yet Another Photon Classification) photon-density-based classification\n",
    "\n",
    "### What is demonstrated\n",
    "\n",
    "* The `icesat2.atl03sp` API is used to perform a SlideRule parallel subsetting request of the Grand Mesa region\n",
    "* The `earthdata.cmr` API's is used to find specific ATL03 granules corresponding to the Grand Mesa region\n",
    "* The `matplotlib` package is used to plot the ATL03 data subset by SlideRule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fb9beaf5-794d-457e-8911-f65158139634",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\") # suppress warnings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54813d4f-f6bd-4fe0-b252-320a1e6804c0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sliderule import sliderule, icesat2, earthdata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "13d96296",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "icesat2.init(\"slideruleearth.io\", verbose=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76b29139-5a28-4744-9436-902d203bc792",
   "metadata": {},
   "source": [
    "## Intro\n",
    "This notebook demonstrates how to use the SlideRule Icesat-2 API to retrieve ATL03 data with two different classifications, one based on the external ATL08-product classifications, designed to distinguish between vegetation and ground returns, and the other based on the experimental YAPC (Yet Another Photon Class) algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c29474c-b649-4739-a3b7-308c1218f76b",
   "metadata": {},
   "source": [
    "### Retrieve ATL03 elevations with ATL08 classifications\n",
    "\n",
    "define a polygon to encompass Grand Mesa, and pick an ATL03 granule that has good coverage over the top of the mesa.  Note that this granule was captured at night, under clear-sky conditions.  Other granules are unlikely to have results as clear s these."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3605d245",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "# build sliderule parameters for ATL03 subsetting request\n",
    "# SRT_LAND = 0\n",
    "# SRT_OCEAN = 1\n",
    "# SRT_SEA_ICE = 2\n",
    "# SRT_LAND_ICE = 3\n",
    "# SRT_INLAND_WATER = 4\n",
    "parms = {\n",
    "    # processing parameters\n",
    "    \"srt\": icesat2.SRT_LAND,\n",
    "    \"len\": 20,\n",
    "    \"res\": 20,\n",
    "    # classification and checks\n",
    "    # still return photon segments that fail checks\n",
    "    \"pass_invalid\": True, \n",
    "    # all photons\n",
    "    \"cnf\": -2, \n",
    "    # all land classification flags\n",
    "    \"atl08_class\": [\"atl08_noise\", \"atl08_ground\", \"atl08_canopy\", \"atl08_top_of_canopy\", \"atl08_unclassified\"],\n",
    "    # all photons\n",
    "    \"yapc\": dict(knn=0, win_h=6, win_x=11, min_ph=4, score=0), \n",
    "}\n",
    "\n",
    "# ICESat-2 data release\n",
    "release = '006'\n",
    "# region of interest\n",
    "poly = [{'lat': 39.34603060272382, 'lon': -108.40601489205419},\n",
    "  {'lat': 39.32770853617356, 'lon': -107.68485163209928},\n",
    "  {'lat': 38.770676045922684, 'lon': -107.71081820956682},\n",
    "  {'lat': 38.788639821001155, 'lon': -108.42635020791396},\n",
    "  {'lat': 39.34603060272382, 'lon': -108.40601489205419}]\n",
    "# time bounds for CMR query\n",
    "time_start = '2019-11-14'\n",
    "time_end = '2019-11-15'\n",
    "\n",
    "# find granule for each region of interest\n",
    "granules_list = earthdata.cmr(short_name='ATL03', polygon=poly, time_start=time_start, time_end=time_end, version=release)\n",
    "\n",
    "# create an empty geodataframe\n",
    "parms[\"poly\"] = poly\n",
    "gdf = icesat2.atl03sp(parms, resources=granules_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e7c90b1-579a-4f4e-a6ab-0c41c1195963",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "gdf.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ac6a168-a8df-46ed-bdcb-64633dedf0da",
   "metadata": {},
   "source": [
    "### Reduce GeoDataFrame to plot a single beam\n",
    "- Convert coordinate reference system to compound projection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38a6f5b5-08f0-45ff-8d8d-ab392d653a4d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def reduce_dataframe(gdf, RGT=None, GT=None, track=None, pair=None, cycle=None, beam='', crs=4326):\n",
    "    # convert coordinate reference system\n",
    "    D3 = gdf.to_crs(crs)\n",
    "    # reduce to reference ground track\n",
    "    if RGT is not None:\n",
    "        D3 = D3[D3[\"rgt\"] == RGT]\n",
    "    # reduce to ground track (gt[123][lr]), track ([123]), or pair (l=0, r=1) \n",
    "    gtlookup = {icesat2.GT1L: 1, icesat2.GT1R: 1, icesat2.GT2L: 2, icesat2.GT2R: 2, icesat2.GT3L: 3, icesat2.GT3R: 3}\n",
    "    pairlookup = {icesat2.GT1L: 0, icesat2.GT1R: 1, icesat2.GT2L: 0, icesat2.GT2R: 1, icesat2.GT3L: 0, icesat2.GT3R: 1}\n",
    "    if GT is not None:\n",
    "        D3 = D3[(D3[\"track\"] == gtlookup[GT]) & (D3[\"pair\"] == pairlookup[GT])]\n",
    "    if track is not None:\n",
    "        D3 = D3[D3[\"track\"] == track]\n",
    "    if pair is not None:\n",
    "        D3 = D3[D3[\"pair\"] == pair]\n",
    "    # reduce to weak or strong beams\n",
    "    # tested on cycle 11, where the strong beam in the pair matches the spacecraft orientation.\n",
    "    # Need to check on other cycles\n",
    "    if (beam == 'strong'):\n",
    "        D3 = D3[D3['sc_orient'] == D3['pair']]\n",
    "    elif (beam == 'weak'):\n",
    "        D3 = D3[D3['sc_orient'] != D3['pair']]\n",
    "    # reduce to cycle\n",
    "    if cycle is not None:\n",
    "        D3 = D3[D3[\"cycle\"] == cycle]\n",
    "    # otherwise, return both beams\n",
    "    return D3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5afec4e1-44c0-44b6-b9dd-7eb28375dfce",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "beam_type = 'strong'\n",
    "project_srs = \"EPSG:26912+EPSG:5703\"\n",
    "D3 = reduce_dataframe(gdf, RGT=737, track=1, beam='strong', crs=project_srs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "92ec1a67-4581-42d1-9f50-8a16586cea40",
   "metadata": {},
   "source": [
    "### Inspect Coordinate Reference System"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9dce7e2a-a9bf-4028-bfcb-628881cddcc5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "D3.crs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15dd2545-818a-4cd3-aa1b-85e90be1e065",
   "metadata": {},
   "source": [
    "### Plot the ATL08 classifications"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80e69cb2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=[8,6])\n",
    "\n",
    "colors={0:['gray', 'noise'], \n",
    "        4:['gray','unclassified'],  \n",
    "        2:['green','canopy'], \n",
    "        3:['lime', 'canopy_top'], \n",
    "        1:['brown', 'ground']}\n",
    "d0=np.min(D3['segment_dist'])\n",
    "for class_val, color_name in colors.items():\n",
    "    ii=D3['atl08_class']==class_val\n",
    "    plt.plot(D3['segment_dist'][ii]+D3['x_atc'][ii]-d0, D3['height'][ii],'o', \n",
    "         markersize=1, color=color_name[0], label=color_name[1])\n",
    "hl=plt.legend(loc=3, frameon=False, markerscale=5)\n",
    "plt.gca().set_xlim([26000, 30000])\n",
    "plt.gca().set_ylim([2950, 3150])\n",
    "\n",
    "plt.ylabel('height, m')\n",
    "plt.xlabel('$x_{ATC}$, m');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a7b578b0-393b-4812-bee7-08749311d4b4",
   "metadata": {},
   "source": [
    "### Plot the YAPC classifications"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ec925a6-161c-4eb8-bfd8-4aa0a785f124",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=[10,6])\n",
    "\n",
    "d0=np.min(D3['segment_dist'])\n",
    "ii=np.argsort(D3['yapc_score'])\n",
    "plt.scatter(D3['segment_dist'][ii]+D3['x_atc'][ii]-d0,\n",
    "    D3['height'][ii],2, c=D3['yapc_score'][ii],\n",
    "    vmin=100, vmax=255, cmap='plasma_r')\n",
    "plt.colorbar(label='YAPC score')\n",
    "plt.gca().set_xlim([26000, 30000])\n",
    "plt.gca().set_ylim([2950, 3150])\n",
    "\n",
    "plt.ylabel('height, m')\n",
    "plt.xlabel('$x_{ATC}$, m');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0ba4dfff-e07a-4008-820a-5a6a24f33e62",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}