{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bd51a330-73c7-494a-b3d9-3f6f07935604",
   "metadata": {},
   "source": [
    "## ArcticDEM Strips Example\n",
    "\n",
    "### Purpose\n",
    "Demonstrate how to work with individual strips when sampling ArcticDEM at ATL06-SR points\n",
    "\n",
    "### Prerequisites\n",
    "1. Access to the PGC S3 bucket `pgc-opendata-dems`\n",
    "2. `gdalinfo` tool installed local to jupyter lab"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a12e5062-ed39-4694-ab2a-53ba804f34ec",
   "metadata": {},
   "source": [
    "#### Import Packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9af774d2-b404-4c52-b932-83992803675f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\") # suppress warnings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7be24a2-d09a-4ed3-8252-e010bd5be5b9",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import sliderule\n",
    "from sliderule import icesat2\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import pyproj\n",
    "import re\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ef5f6de-1454-4c50-82f9-d01c252aede9",
   "metadata": {},
   "source": [
    "#### Initialize Python Client"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48c7b8ca-74bf-4998-965d-d8b871a10c51",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "icesat2.init(\"slideruleearth.io\", verbose=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "751c7ba4-3985-4d84-a98b-89eae6346dd7",
   "metadata": {},
   "source": [
    "#### Build Region of Interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04333ea1-3564-4657-a392-bee8d0c5de62",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "xy0=np.array([  -73000., -2683000.])\n",
    "transformer = pyproj.Transformer.from_crs(3413, 4326)\n",
    "xyB=[xy0[0]+np.array([-1, 1, 1, -1, -1])*1.e4, xy0[1]+np.array([-1, -1, 1, 1, -1])*1.e4]\n",
    "llB=transformer.transform(*xyB)\n",
    "poly=[{'lat':lat,'lon':lon} for lat, lon in zip(*llB)]\n",
    "plist = []\n",
    "for p in poly:\n",
    "    plist += p[\"lat\"], \n",
    "    plist += p[\"lon\"],\n",
    "region_of_interest = sliderule.toregion(plist)\n",
    "region_of_interest[\"gdf\"].plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea351ded-8069-4f0c-a4bc-c41854e5f583",
   "metadata": {},
   "source": [
    "#### Make Processing Request\n",
    "ATL06-SR request includes the `samples` parameter to specify that ArcticDEM Strips dataset should be sampled at each generated ATL06 elevation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79ad73d6-ece4-4ebd-ac90-77c4d7cf006f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "parms = { \"poly\": region_of_interest[\"poly\"],\n",
    "          \"cnf\": \"atl03_high\",\n",
    "          \"ats\": 10.0,\n",
    "          \"cnt\": 5,\n",
    "          \"len\": 40.0,\n",
    "          \"res\": 120.0,\n",
    "          \"rgt\": 658,\n",
    "          \"time_start\":'2020-01-01',\n",
    "          \"time_end\":'2021-01-01',\n",
    "          \"samples\": {\"strips\": {\"asset\": \"arcticdem-strips\", \"with_flags\": True}} }\n",
    "gdf = icesat2.atl06p(parms)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cf78cf33-9ed1-402f-b2c9-be8604ed823e",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### Print Out File Directory\n",
    "When a GeoDataFrame includes samples from rasters, each sample value has a file id that is used to look up the file name of the source raster for that value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e330379c-3890-4314-bb74-9780e549d724",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "gdf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19c9060d-c502-425c-99c4-7d55a9780ef5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "gdf.attrs['file_directory']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6dca9849-713f-4c7a-a68e-e6b782fac3b0",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### Pull Out Bounding Box of Raster\n",
    "This step requires AWS credentials to be able to access S3 and `gdalinfo` be installed on the host machine to read the bounding box for each raster sampled by SlideRule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a97e1d1-ff1c-4068-a305-04959faf5d49",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# helper functions\n",
    "def getXY(line):\n",
    "    line = line.replace(\"(\",\"$\")\n",
    "    line = line.replace(\")\",\"$\")\n",
    "    point = line.split(\"$\")[1]\n",
    "    coord = point.split(\",\")\n",
    "    x = float(coord[0].strip())\n",
    "    y = float(coord[1].strip())\n",
    "    return x, y\n",
    "\n",
    "def getLonLat(line):\n",
    "    line = line.replace(\"(\",\"$\")\n",
    "    line = line.replace(\")\",\"$\")\n",
    "    point = line.split(\"$\")[3]\n",
    "    coord = point.split(\",\")\n",
    "    deg, minutes, seconds, direction = re.split('[d\\'\"]', coord[1].strip())\n",
    "    lon = (float(deg) + float(minutes)/60 + float(seconds)/(60*60)) * (-1 if direction in ['W', 'S'] else 1)\n",
    "    deg, minutes, seconds, direction = re.split('[d\\'\"]', coord[0].strip())\n",
    "    lat = (float(deg) + float(minutes)/60 + float(seconds)/(60*60)) * (-1 if direction in ['W', 'S'] else 1)\n",
    "    return [lon, lat]\n",
    "\n",
    "def getBB(dem):\n",
    "    os.system(\"gdalinfo {} > /tmp/r.txt\".format(dem))\n",
    "    with open(\"/tmp/r.txt\", \"r\") as file:\n",
    "        lines = file.readlines()\n",
    "        for line in lines:\n",
    "            if \"Upper Left\" in line:\n",
    "                ul = getLonLat(line)\n",
    "            elif \"Lower Left\" in line:\n",
    "                ll = getLonLat(line)\n",
    "            elif \"Upper Right\" in line:\n",
    "                ur = getLonLat(line)\n",
    "            elif \"Lower Right\" in line:\n",
    "                lr = getLonLat(line)\n",
    "    return ul + ll + lr + ur + ul"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce2fff54-c837-48d1-96c1-a396f01b54ff",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "file_dict = {}\n",
    "for file_id, file_name in gdf.attrs['file_directory'].items():\n",
    "    if \"bitmask\" in file_name:\n",
    "        continue\n",
    "    if file_name not in file_dict:\n",
    "        file_dict[file_name] = {\"ids\": []}\n",
    "    file_dict[file_name][\"ids\"].append(file_id)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8fd83146-6f4c-4ee1-91b6-6187ea2f5646",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# get boundaries for each raster\n",
    "for file_name in file_dict:\n",
    "    print(\"Retrieving raster info for:\", file_name)\n",
    "    rlist = getBB(file_name)\n",
    "    file_dict[file_name][\"region\"] = sliderule.toregion(rlist)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b65bb3d5-de72-4d03-aada-ffd394741f6a",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### Pull Out Individual DEM Values and Put in Separate Columns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "478f0a19-2aa7-44ea-8c81-3b4b0c8657da",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def getValue(x, file_ids):\n",
    "    for file_id in file_ids:\n",
    "        l = np.where(x['strips.file_id'] == file_id)[0]\n",
    "        if len(l) == 1:\n",
    "            return x['strips.value'][l[0]]\n",
    "    return None\n",
    "sampled_data = gdf[gdf['strips.time'].notnull()]\n",
    "for file_name in file_dict:\n",
    "    sampled_data[file_name] = sampled_data.apply(lambda x: getValue(x, file_dict[file_name][\"ids\"]), axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "574fb9af-6db6-44e2-9a7a-1fab84a15f51",
   "metadata": {},
   "source": [
    "#### Plot Overlays of Boundaries and Returns"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec498d91-77dd-415b-b0fe-8125d0bd6655",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "with warnings.catch_warnings():\n",
    "    warnings.simplefilter(\"ignore\")\n",
    "    fig = plt.figure(num=None, figsize=(24, 24))\n",
    "    region_lons = [p[\"lon\"] for p in region_of_interest[\"poly\"]]\n",
    "    region_lats = [p[\"lat\"] for p in region_of_interest[\"poly\"]]\n",
    "    ax = {}\n",
    "    k = 0\n",
    "    for file_name in file_dict:\n",
    "        raster = file_dict[file_name]\n",
    "        raster_lons = [p[\"lon\"] for p in raster[\"region\"][\"poly\"]]\n",
    "        raster_lats = [p[\"lat\"] for p in raster[\"region\"][\"poly\"]]\n",
    "        plot_data = sampled_data[sampled_data[file_name].notnull()]\n",
    "        ax[k] = plt.subplot(5,4,k+1)\n",
    "        gdf.plot(ax=ax[k], column='h_mean', color='y', markersize=0.5)\n",
    "        plot_data.plot(ax=ax[k], column='h_mean', color='b', markersize=0.5)\n",
    "        ax[k].plot(region_lons, region_lats, linewidth=1.5, color='r', zorder=2)\n",
    "        ax[k].plot(raster_lons, raster_lats, linewidth=1.5, color='g', zorder=2)\n",
    "        k += 1\n",
    "    plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1addca2-fdbb-4e37-ae3e-417652e1f119",
   "metadata": {},
   "source": [
    "#### Plot the Different ArcticDEM Values against the SlideRule ATL06-SR Values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f98551c9-c76b-4469-a6df-acf2c2c375ef",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Select DEM File ID\n",
    "file_name = list(file_dict.keys())[0]\n",
    "\n",
    "# Setup Plot\n",
    "fig,ax = plt.subplots(num=None, figsize=(10, 8))\n",
    "ax.set_title(\"SlideRule vs. ArcticDEM Elevations\")\n",
    "ax.set_xlabel('distance (m)')\n",
    "ax.set_ylabel('height (m)')\n",
    "legend_elements = []\n",
    "\n",
    "# Filter Data to Plot\n",
    "plot_data = sampled_data[sampled_data[file_name].notnull()]\n",
    "\n",
    "# Set X Axis\n",
    "x_axis = plot_data[\"x_atc\"]\n",
    "\n",
    "# Plot SlideRule ATL06 Elevations\n",
    "sc1 = ax.scatter(x_axis, plot_data[\"h_mean\"].values, c='red', s=2.5)\n",
    "legend_elements.append(matplotlib.lines.Line2D([0], [0], color='red', lw=6, label='ATL06-SR'))\n",
    "\n",
    "# Plot ArcticDEM Elevations\n",
    "sc2 = ax.scatter(x_axis, plot_data[file_name].values, c='blue', s=2.5)\n",
    "legend_elements.append(matplotlib.lines.Line2D([0], [0], color='blue', lw=6, label='ArcticDEM'))\n",
    "\n",
    "# Display Legend\n",
    "lgd = ax.legend(handles=legend_elements, loc=3, frameon=True)\n",
    "lgd.get_frame().set_alpha(1.0)\n",
    "lgd.get_frame().set_edgecolor('white')\n",
    "\n",
    "# Show Plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5270acb7-a1f8-4fce-bd16-91c2c055816a",
   "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
}