{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ddce2704-937f-46be-a0c8-fd402cc50037",
   "metadata": {},
   "source": [
    "# Generating ATL03 photon classifications using ATL08 and YAPC\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.atl03x` 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": "e94bd36c",
   "metadata": {},
   "outputs": [],
   "source": [
    "sliderule.init(verbose=True)"
   ]
  },
  {
   "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",
    "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",
    "\n",
    "# region of interest\n",
    "poly = [\n",
    "  {'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",
    "]\n",
    "\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 geodataframe\n",
    "gdf = sliderule.run(\"atl03x\", parms, aoi=poly, resources=granules_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e7c90b1-579a-4f4e-a6ab-0c41c1195963",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "gdf"
   ]
  },
  {
   "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": "7d0bd1dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38a6f5b5-08f0-45ff-8d8d-ab392d653a4d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def reduce_dataframe(gdf, RGT=None, GT=None, spot=None, cycle=None, crs=4326):\n",
    "    D3 = gdf.to_crs(crs) # convert coordinate reference system\n",
    "    if RGT is not None:\n",
    "        D3 = D3[D3[\"rgt\"] == RGT]\n",
    "    if GT is not None:\n",
    "        D3 = D3[D3[\"gt\"] == GT]\n",
    "    if spot is not None:\n",
    "        D3 = D3[D3[\"spot\"] == spot]\n",
    "    if cycle is not None:\n",
    "        D3 = D3[D3[\"cycle\"] == cycle]\n",
    "    return D3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5afec4e1-44c0-44b6-b9dd-7eb28375dfce",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "project_srs = \"EPSG:26912+EPSG:5703\"\n",
    "D3 = reduce_dataframe(gdf, RGT=737, spot=1, 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['x_atc'])\n",
    "for class_val, color_name in colors.items():\n",
    "    ii=D3['atl08_class']==class_val\n",
    "    plt.plot(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([25000, 30000])\n",
    "plt.gca().set_ylim([3050, 3300])\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['x_atc'])\n",
    "ii=np.argsort(D3['yapc_score'])\n",
    "plt.scatter(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([25000, 30000])\n",
    "plt.gca().set_ylim([3050, 3300])\n",
    "\n",
    "plt.ylabel('height, m')\n",
    "plt.xlabel('$x_{ATC}$, m');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "753b9de0",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "sliderule_notebook",
   "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.14.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}