SDAP_visualizations.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ae58a5d3-82b9-428a-875c-23dcbd20b7f7",
   "metadata": {},
   "source": [
    "### The STV-FIS effort is funded by the NASA Decadal Survey Incubation program - **NNH21ZDA001N-DSI**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e523f06-fb9c-45d2-bf17-5c9abacfbc23",
   "metadata": {},
   "source": [
    "# SDAP Tomogram & LIDAR Subsetting & Visualization Endpoints"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4945ed7d-a5e1-43eb-8258-c8403aee031d",
   "metadata": {},
   "source": [
    "This notebook showcases some vizualization capabilities built directly into SDAP for the purposes of visuzlizing SAR, LIDAR, and ML-Fusion data. The following demo products were used:\n",
    "- SAR: Work in progress geolocated product (not publicly availabe). Kings Canyon, CA; Lope Nat'l Park, Gabon (HH, HV, VV); Rabi Forest, Gabon (HH, HV, VV); Redwood Nat'l Park, CA (HH, HV, VV | L- & P-band)\n",
    "- LIDAR: [This dataset](https://daac.ornl.gov/ABOVE/guides/ABoVE_LVIS_VegetationStructure.html) (ABoVE: LVIS L3 Gridded Vegetation Structure across North America, 2017 and 2019)\n",
    "- Kauri: An ML fusion product using SAR and LIDAR data to predict RH98 (canopy height)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9172578d-d745-4ff2-a031-850c64e6738b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "import os\n",
    "import sys\n",
    "from base64 import b64encode\n",
    "from copy import deepcopy\n",
    "from datetime import datetime\n",
    "from io import BytesIO, StringIO\n",
    "from math import floor, ceil\n",
    "from zipfile import ZipFile\n",
    "\n",
    "import fiona\n",
    "import ipywidgets as widgets\n",
    "import matplotlib.pyplot as plt\n",
    "import geopandas as gpd\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import requests\n",
    "import urllib3\n",
    "import xarray as xr\n",
    "from ipyleaflet import DrawControl, Map, basemaps, GeoData, ImageOverlay\n",
    "from matplotlib.cm import ScalarMappable\n",
    "from matplotlib.colors import Normalize\n",
    "from requests.compat import urlencode\n",
    "from shapely import from_geojson\n",
    "from shapely.geometry import Point\n",
    "from IPython.display import display\n",
    "from IPython.display import Image as IpImage\n",
    "from PIL import Image\n",
    "\n",
    "fiona.drvsupport.supported_drivers[\"LIBKML\"] = \"rw\"\n",
    "urllib3.disable_warnings(urllib3.exceptions.InsecureRequestWarning)\n",
    "\n",
    "# NOTE: Disabling verification due to certificate issues on my local machine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "33a4daf5-2c03-4a8c-80cf-a3acb52c07ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "BASE_URL = 'https://ideas-digitaltwin.jpl.nasa.gov'\n",
    "NEXUS_URL = f'{BASE_URL}/nexus'\n",
    "\n",
    "DATASETS_URL = f'{NEXUS_URL}/list'\n",
    "\n",
    "TOMOGRAM_LAT_URL = f'{NEXUS_URL}/tomogram/latitude'\n",
    "TOMOGRAM_LON_URL = f'{NEXUS_URL}/tomogram/longitude'\n",
    "TOMOGRAM_CUS_URL = f'{NEXUS_URL}/tomogram/profile'\n",
    "TOMOGRAM_3D_URL = f'{NEXUS_URL}/tomogram/3d'\n",
    "\n",
    "LIDAR_URL = f'{NEXUS_URL}/stv/lidar'\n",
    "\n",
    "TOMOGRAM_DS_KC = 'Kings_Canyon_Tomogram'\n",
    "TOMOGRAM_DS_Lope = 'Lope_Tomogram_{}'\n",
    "TOMOGRAM_DS_Rabi = 'Rabi_Tomogram_{}'\n",
    "TOMOGRAM_DS_Redwood_Lband = 'Redwood_Tomogram_L_{}'\n",
    "TOMOGRAM_DS_Redwood_Pband = 'Redwood_Tomogram_P_{}'\n",
    "\n",
    "LIDAR_DS = 'LVIS_ABoVE'\n",
    "\n",
    "KAURI_SAMPLE_DS = 'Kauri_sample_output'\n",
    "\n",
    "TOMO_MAP_DATASETS = { \n",
    "    TOMOGRAM_DS_Rabi.format('HH'): (\n",
    "        [\n",
    "            'Gabon_LVIS_DSM', \n",
    "            'Rabi_HybridSAR_L_H5_AfriSAR_DSM', \n",
    "            'rabi_L_KauriCNN_canopy_elevation_5x5_200epochs', \n",
    "        ], \n",
    "        [\n",
    "            'Gabon_LVIS_ZG', \n",
    "            'rabi_L_KauriCNN_ground_elevation_dem_5x5_60epochs', \n",
    "        ]\n",
    "    ),\n",
    "    TOMOGRAM_DS_Rabi.format('HV'): (\n",
    "        [\n",
    "            'Gabon_LVIS_DSM', \n",
    "            'Rabi_HybridSAR_L_H5_AfriSAR_DSM', \n",
    "            'rabi_L_KauriCNN_canopy_elevation_5x5_200epochs', \n",
    "        ], \n",
    "        [\n",
    "            'Gabon_LVIS_ZG', \n",
    "            'rabi_L_KauriCNN_ground_elevation_dem_5x5_60epochs', \n",
    "        ]\n",
    "    ),\n",
    "    TOMOGRAM_DS_Rabi.format('VV'): (\n",
    "        [\n",
    "            'Gabon_LVIS_DSM', \n",
    "            'Rabi_HybridSAR_L_H5_AfriSAR_DSM', \n",
    "            'rabi_L_KauriCNN_canopy_elevation_5x5_200epochs', \n",
    "        ], \n",
    "        [\n",
    "            'Gabon_LVIS_ZG', \n",
    "            'rabi_L_KauriCNN_ground_elevation_dem_5x5_60epochs', \n",
    "        ]\n",
    "    ),\n",
    "    TOMOGRAM_DS_Lope.format('HH'): (\n",
    "        ['Gabon_LVIS_DSM', 'Lope_HybridSAR_L_H5_AfriSAR_DSM'], \n",
    "        ['Gabon_LVIS_ZG', 'Lope_Reference_DTM']\n",
    "    ),\n",
    "    TOMOGRAM_DS_Lope.format('HV'): (\n",
    "        ['Gabon_LVIS_DSM', 'Lope_HybridSAR_L_H5_AfriSAR_DSM'], \n",
    "        ['Gabon_LVIS_ZG', 'Lope_Reference_DTM']\n",
    "    ),\n",
    "    TOMOGRAM_DS_Lope.format('VV'): (\n",
    "        ['Gabon_LVIS_DSM', 'Lope_HybridSAR_L_H5_AfriSAR_DSM'], \n",
    "        ['Gabon_LVIS_ZG', 'Lope_Reference_DTM']\n",
    "    ),\n",
    "    (TOMOGRAM_DS_Rabi + '_NO_DEM').format('HH'): (\n",
    "        [\n",
    "            'Gabon_LVIS_RH098', \n",
    "            'Rabi_HybridSAR_L_H5_AfriSAR_CHM', \n",
    "            'rabi_L_KauriCNN_canopy_height_5x5_200epochs'\n",
    "        ], \n",
    "        [\n",
    "            'rabi_L_KauriCNN_ground_elevation_bias_5x5_60epochs'\n",
    "        ]),\n",
    "    (TOMOGRAM_DS_Rabi + '_NO_DEM').format('HV'): (\n",
    "        [\n",
    "            'Gabon_LVIS_RH098', \n",
    "            'Rabi_HybridSAR_L_H5_AfriSAR_CHM', \n",
    "            'rabi_L_KauriCNN_canopy_height_5x5_200epochs'\n",
    "        ], \n",
    "        [\n",
    "            'rabi_L_KauriCNN_ground_elevation_bias_5x5_60epochs']),\n",
    "    (TOMOGRAM_DS_Rabi + '_NO_DEM').format('VV'): (\n",
    "        [\n",
    "            'Gabon_LVIS_RH098', \n",
    "            'Rabi_HybridSAR_L_H5_AfriSAR_CHM', \n",
    "            'rabi_L_KauriCNN_canopy_height_5x5_200epochs'\n",
    "        ], \n",
    "        [\n",
    "            'rabi_L_KauriCNN_ground_elevation_bias_5x5_60epochs'\n",
    "        ]),\n",
    "    (TOMOGRAM_DS_Lope + '_NO_DEM').format('HH'): (\n",
    "        [\n",
    "            'Gabon_LVIS_RH098', \n",
    "            'Lope_HybridSAR_L_H5_AfriSAR_CHM' \n",
    "        ], \n",
    "        []\n",
    "    ),\n",
    "    (TOMOGRAM_DS_Lope + '_NO_DEM').format('HV'): (\n",
    "        [\n",
    "            'Gabon_LVIS_RH098', \n",
    "            'Lope_HybridSAR_L_H5_AfriSAR_CHM'\n",
    "        ], \n",
    "        []\n",
    "    ),\n",
    "    (TOMOGRAM_DS_Lope + '_NO_DEM').format('VV'): (\n",
    "        [\n",
    "            'Gabon_LVIS_RH098', \n",
    "            'Lope_HybridSAR_L_H5_AfriSAR_CHM'\n",
    "        ], \n",
    "        []\n",
    "    ),\n",
    "}\n",
    "\n",
    "\n",
    "def get(url, params=None):\n",
    "    t = datetime.now()\n",
    "        \n",
    "    print(f'{url}', end='')\n",
    "\n",
    "    if params is not None:\n",
    "        print(f'?{urlencode(params)}')\n",
    "    else:\n",
    "        print()\n",
    "        params = dict()\n",
    "    \n",
    "    print('Waiting for response from SDAP...', end='')\n",
    "    response = requests.get(url, params=params, verify=False)\n",
    "    td = datetime.now() - t\n",
    "\n",
    "    if response.status_code == 200:\n",
    "        print(f' Done in {td}')\n",
    "        return response\n",
    "    elif response.status_code == 200:\n",
    "        print(f' Timed out after {td}')\n",
    "        return None\n",
    "    else:\n",
    "        print(f' Failed in {td}')\n",
    "        print(response.text)\n",
    "        response.raise_for_status()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7e467b15-f7b6-4693-a26c-0dfdbb701ae5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found 39 datasets described in this notebook\n",
      "https://ideas-digitaltwin.jpl.nasa.gov/nexus/list\n",
      "Waiting for response from SDAP... Done in 0:00:00.200912\n",
      "Found 334 datasets available in https://ideas-digitaltwin.jpl.nasa.gov\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: Found 19 datasets used in this notebook that are missing from https://ideas-digitaltwin.jpl.nasa.gov\n",
      "\t- LVIS_ABoVE\n",
      "\t- Lope_Tomogram_HH\n",
      "\t- Lope_Tomogram_HH_NO_DEM\n",
      "\t- Lope_Tomogram_HV\n",
      "\t- Lope_Tomogram_HV_NO_DEM\n",
      "\t- Lope_Tomogram_VV\n",
      "\t- Lope_Tomogram_VV_NO_DEM\n",
      "\t- Redwood_Tomogram_L_HH\n",
      "\t- Redwood_Tomogram_L_HH_NO_DEM\n",
      "\t- Redwood_Tomogram_L_HV\n",
      "\t- Redwood_Tomogram_L_HV_NO_DEM\n",
      "\t- Redwood_Tomogram_L_VV\n",
      "\t- Redwood_Tomogram_L_VV_NO_DEM\n",
      "\t- Redwood_Tomogram_P_HH\n",
      "\t- Redwood_Tomogram_P_HH_NO_DEM\n",
      "\t- Redwood_Tomogram_P_HV\n",
      "\t- Redwood_Tomogram_P_HV_NO_DEM\n",
      "\t- Redwood_Tomogram_P_VV\n",
      "\t- Redwood_Tomogram_P_VV_NO_DEM\n",
      "\n",
      "Make note of these and avoid using queries that reference them\n"
     ]
    }
   ],
   "source": [
    "# Audit available datasets\n",
    "\n",
    "datasets = []\n",
    "tomo_datasets = [TOMOGRAM_DS_KC]\n",
    "\n",
    "for tomo_site in [TOMOGRAM_DS_Lope, TOMOGRAM_DS_Rabi, TOMOGRAM_DS_Redwood_Lband, TOMOGRAM_DS_Redwood_Pband]:\n",
    "    for pol in ['HH', 'HV', 'VV']:\n",
    "        tomo_datasets.extend([tomo_site.format(pol), tomo_site.format(pol) + \"_NO_DEM\"])\n",
    "\n",
    "datasets.extend(tomo_datasets)\n",
    "\n",
    "datasets.append(LIDAR_DS)\n",
    "\n",
    "maps = []\n",
    "\n",
    "for v in TOMO_MAP_DATASETS.values():\n",
    "    for p in v:\n",
    "        maps.extend(p)\n",
    "\n",
    "datasets.extend(set(maps))\n",
    "\n",
    "datasets.append(KAURI_SAMPLE_DS)\n",
    "\n",
    "datasets.sort()\n",
    "\n",
    "print(f'Found {len(datasets)} datasets described in this notebook')\n",
    "\n",
    "available_datasets = [ds['shortName'] for ds in get(DATASETS_URL).json()]\n",
    "\n",
    "print(f'Found {len(available_datasets)} datasets available in {BASE_URL}')\n",
    "\n",
    "missing = list(set(datasets).difference(set(available_datasets)))\n",
    "\n",
    "if len(missing) > 0:\n",
    "    print(f'WARNING: Found {len(missing)} datasets used in this notebook that are missing from {BASE_URL}', file=sys.stderr)\n",
    "    for m in sorted(missing):\n",
    "        print(f'\\t- {m}', file=sys.stderr)\n",
    "    print(file=sys.stderr)\n",
    "    print('Make note of these and avoid using queries that reference them', file=sys.stderr)\n",
    "else:\n",
    "    print('All datasets are present')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf7a5d70-1832-443e-b35e-b0f01f112382",
   "metadata": {},
   "source": [
    "## Tomograms"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57a61452-a64b-4f18-b796-fe0c46edae43",
   "metadata": {},
   "source": [
    "### Tomogram Profile Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "beffa2e8-4a55-4c6c-8f2c-8a196d82fe45",
   "metadata": {},
   "outputs": [],
   "source": [
    "sites = [\n",
    "    ('Kings Canyon National Park, US', 'kcnp'),\n",
    "    ('Lope National Park, GA', 'lope'),\n",
    "    ('Rabi Forest, GA', 'rabi'),\n",
    "    ('Redwood National Park, US (L-band)', 'redl'),\n",
    "    ('Redwood National Park, US (P-band)', 'redp'),\n",
    "]\n",
    "\n",
    "site_sel = widgets.Dropdown(options=sites, description='Site: ', disabled=False, layout = widgets.Layout(width='350px'))\n",
    "pol_sel = widgets.Dropdown(options=['HH', 'HV', 'VV'], description='Pol: ', disabled=False, layout = widgets.Layout(width='350px'))\n",
    "dim_sel = widgets.Dropdown(options=['Latitude', 'Longitude', 'Custom'], description='Direction: ', disabled=False, layout = widgets.Layout(width='350px'))\n",
    "dem_sel = widgets.Checkbox(value=True, disabled=False, description='Add DEM')\n",
    "peaks = widgets.Checkbox(value=False, disabled=False, description='Show peaks in plot')\n",
    "\n",
    "\n",
    "display(site_sel)\n",
    "display(pol_sel)\n",
    "display(dim_sel)\n",
    "display(dem_sel)\n",
    "display(peaks)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "abe3904c-6ed5-4306-9582-10af98d23a37",
   "metadata": {},
   "outputs": [],
   "source": [
    "PARAMS = {\n",
    "    'kcnp': dict(ds=TOMOGRAM_DS_KC,            bounds=dict(min_lon=-119.60,  max_lon=-118.75,  min_lat=36.98,  max_lat=37.15,  min_elev=66.84,  max_elev=3732.20, min_elev_nd=-60,  max_elev_nd=60), kmzpng='kmz/kingsF_09433_22046_021_220929_PL09043020_05_CX_01'),\n",
    "    'lope': dict(ds=TOMOGRAM_DS_Lope,          bounds=dict(min_lon=11.43,    max_lon=11.83,    min_lat=-0.34,  max_lat=0.07,   min_elev=2.5,    max_elev=717.87,  min_elev_nd=-102, max_elev_nd=56), kmzpng='kmz/lopenp_14043_16008_009_160225_L090_CX_129A_02'),\n",
    "    'rabi': dict(ds=TOMOGRAM_DS_Rabi,          bounds=dict(min_lon=9.70,     max_lon=10.11,    min_lat=-2.05,  max_lat=-1.65,  min_elev=-94.90, max_elev=236.65,  min_elev_nd=-102, max_elev_nd=56), kmzpng='kmz/rabifo_13049_16010_010_160228_L090_CX_129A_03'),\n",
    "    'redl': dict(ds=TOMOGRAM_DS_Redwood_Lband, bounds=dict(min_lon=-124.247, max_lon=-123.903, min_lat=41.571, max_lat=41.888, min_elev=-89.13, max_elev=1147.06, min_elev_nd=-60,  max_elev_nd=60), kmzpng='kmz/redwoo_16912_22043_019_220913_L090_CX_01'),\n",
    "    'redp': dict(ds=TOMOGRAM_DS_Redwood_Pband, bounds=dict(min_lon=-124.247, max_lon=-123.896, min_lat=41.528, max_lat=41.888, min_elev=-89.22, max_elev=1146.80, min_elev_nd=-60,  max_elev_nd=60), kmzpng='kmz/redwoo_16912_22043_019_220913_L090_CX_01'),\n",
    "}\n",
    "\n",
    "params = PARAMS[site_sel.value]\n",
    "\n",
    "if site_sel.value != 'kcnp':\n",
    "    params['ds'] = params['ds'].format(pol_sel.value)\n",
    "    \n",
    "if not dem_sel.value:\n",
    "    params['ds'] = params['ds'] + '_NO_DEM'\n",
    "    params['bounds']['min_elev'] = params['bounds']['min_elev_nd']\n",
    "    params['bounds']['max_elev'] = params['bounds']['max_elev_nd']\n",
    "\n",
    "del params['bounds']['min_elev_nd']\n",
    "del params['bounds']['max_elev_nd']\n",
    "\n",
    "line = {\n",
    "    'start-point': 'invalid',\n",
    "    'end-point': 'invalid'\n",
    "}\n",
    "\n",
    "if dim_sel.value in ['Longitude', 'Latitude']:\n",
    "    s = widgets.FloatSlider(\n",
    "        value=((params['bounds']['min_lon'] + params['bounds']['max_lon']) / 2) if dim_sel.value == 'Longitude' else ((params['bounds']['min_lat'] + params['bounds']['max_lat']) / 2),\n",
    "        min=params['bounds']['min_lon'] if dim_sel.value == 'Longitude' else params['bounds']['min_lat'],\n",
    "        max=params['bounds']['max_lon'] if dim_sel.value == 'Longitude' else params['bounds']['max_lat'],\n",
    "        step=0.001,\n",
    "        description=f'{dim_sel.value}: ',\n",
    "        disabled=False,\n",
    "        continuous_update=False,\n",
    "        orientation='horizontal',\n",
    "        readout=True,\n",
    "        readout_format='.3f',\n",
    "        layout = widgets.Layout(width='500px')\n",
    "    )\n",
    "    \n",
    "    d = widgets.FloatRangeSlider(\n",
    "        value=[\n",
    "            params['bounds']['min_lat'] if dim_sel.value == 'Longitude' else params['bounds']['min_lon'],\n",
    "            params['bounds']['max_lat'] if dim_sel.value == 'Longitude' else params['bounds']['max_lon']\n",
    "        ],\n",
    "        min=params['bounds']['min_lat'] if dim_sel.value == 'Longitude' else params['bounds']['min_lon'],\n",
    "        max=params['bounds']['max_lat'] if dim_sel.value == 'Longitude' else params['bounds']['max_lon'],\n",
    "        step=0.01,\n",
    "        description='Latitude: ' if dim_sel.value == 'Longitude' else 'Longitude: ',\n",
    "        disabled=False,\n",
    "        continuous_update=False,\n",
    "        orientation='horizontal',\n",
    "        readout=True,\n",
    "        readout_format='.2f',\n",
    "        layout = widgets.Layout(width='500px')\n",
    "    )\n",
    "    \n",
    "    display(s)\n",
    "    display(d)\n",
    "else:\n",
    "    kmz_png = params['kmzpng']\n",
    "\n",
    "    gdf = gpd.read_file(kmz_png + '.kmz')\n",
    "    center = gdf.to_crs('EPSG:4326').dissolve().centroid.get_coordinates().stack()[::-1].tolist()\n",
    "\n",
    "    m = Map(center = center, zoom = 10, min_zoom = 9, max_zoom = 20, basemap=basemaps.Esri.WorldImagery)\n",
    "\n",
    "    draw_control = DrawControl()\n",
    "\n",
    "    draw_control.circlemarker = {}\n",
    "    draw_control.polygon = {}\n",
    "    \n",
    "    line = {\n",
    "        'start-point': 'invalid',\n",
    "        'end-point': 'invalid'\n",
    "    }\n",
    "    \n",
    "    def handle_draw(self, action, geo_json):\n",
    "        try:\n",
    "            geom = from_geojson(json.dumps(geo_json))\n",
    "            points = list(geom.coords)\n",
    "        \n",
    "            if len(points) == 2:\n",
    "                line['start-point'] = Point(points[0])    \n",
    "                line['end-point'] = Point(points[1])\n",
    "            else:\n",
    "                line['start-point'] = 'invalid'\n",
    "                line['end-point'] = 'invalid'\n",
    "        except Exception as err:\n",
    "            line['start-point'] = 'invalid'\n",
    "            line['end-point'] = 'invalid'\n",
    "            line['err'] = repr(err)\n",
    "    \n",
    "    draw_control.on_draw(handle_draw)\n",
    "    \n",
    "    m.add_control(draw_control)\n",
    "    \n",
    "    gdata = GeoData(geo_dataframe=gdf, style={'fillOpacity': 0})\n",
    "    \n",
    "    bounds = gdf.dissolve().to_crs('EPSG:4326').bounds.iloc[0]\n",
    "    bounds = ((bounds.miny, bounds.minx), (bounds.maxy, bounds.maxx))\n",
    "    \n",
    "    image = Image.open(kmz_png + '.png')\n",
    "    buf = BytesIO()\n",
    "    \n",
    "    ext = 'png'\n",
    "    \n",
    "    image.save(buf, ext)\n",
    "    \n",
    "    img_data = b64encode(buf.getvalue())\n",
    "    img_data = img_data.decode('ascii')\n",
    "    url = 'data:image{};base64,'.format(ext) + img_data\n",
    "    \n",
    "    img = ImageOverlay(url=url, bounds=bounds)\n",
    "    \n",
    "    m.add(gdata)\n",
    "    m.add(img)\n",
    "    \n",
    "    display(m)\n",
    "\n",
    "e = widgets.IntRangeSlider(\n",
    "    value=[params['bounds']['min_elev'], params['bounds']['max_elev']],\n",
    "    min=floor(params['bounds']['min_elev']),\n",
    "    max=ceil(params['bounds']['max_elev']),\n",
    "    step=1,\n",
    "    description='Elevation: ',\n",
    "    disabled=False,\n",
    "    continuous_update=False,\n",
    "    orientation='horizontal',\n",
    "    readout=True,\n",
    "    readout_format='d',\n",
    "    layout = widgets.Layout(width='500px')\n",
    ")\n",
    "\n",
    "display(e)\n",
    "    \n",
    "options_c, options_g = deepcopy(TOMO_MAP_DATASETS).get(params['ds'], ([], []))\n",
    "\n",
    "if '--' not in options_c:\n",
    "    options_c.append('--')\n",
    "\n",
    "if '--' not in options_g:\n",
    "    options_g.append('--')\n",
    "\n",
    "c_sel = widgets.Dropdown(options=options_c, description='Primary Canopy height map: ', value='--', disabled=False, layout = widgets.Layout(width='500px'), style={'description_width': 'initial'})\n",
    "g_sel = widgets.Dropdown(options=options_g, description='Primary Ground height map: ', value='--', disabled=False, layout = widgets.Layout(width='500px'), style={'description_width': 'initial'})\n",
    "\n",
    "c2_sel = widgets.Dropdown(options=options_c, description='Secondary Canopy height map: ', value='--', disabled=False, layout = widgets.Layout(width='500px'), style={'description_width': 'initial'})\n",
    "g2_sel = widgets.Dropdown(options=options_g, description='Secondary Ground height map: ', value='--', disabled=False, layout = widgets.Layout(width='500px'), style={'description_width': 'initial'})\n",
    "\n",
    "display(c_sel)\n",
    "display(c2_sel)\n",
    "display(g_sel)\n",
    "display(g2_sel)\n",
    "\n",
    "percentile = widgets.Checkbox(value=False, disabled=False, description='Render by cumulative vertical return power', style={'description_width': 'initial'})\n",
    "display(percentile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "205073ea-de90-4ad6-8107-086ed0380e96",
   "metadata": {},
   "outputs": [],
   "source": [
    "if dim_sel.value in ['Longitude', 'Latitude']:\n",
    "    if dim_sel.value == 'Longitude':\n",
    "        params['longitude'] = s.value\n",
    "        (params['minLat'], params['maxLat']) = d.value\n",
    "        url = TOMOGRAM_LON_URL\n",
    "    else:\n",
    "        params['latitude'] = s.value\n",
    "        (params['minLon'], params['maxLon']) = d.value\n",
    "        url = TOMOGRAM_LAT_URL\n",
    "else:\n",
    "    draw_control.clear_polylines()\n",
    "    if not isinstance(line['start-point'], Point):\n",
    "        raise ValueError('Invalid line selection')\n",
    "\n",
    "    start = line['start-point']\n",
    "    end = line['end-point']\n",
    "    \n",
    "    params['p1'] = f'{start.x},{start.y}'\n",
    "    params['p2'] = f'{end.x},{end.y}'\n",
    "    \n",
    "    url = TOMOGRAM_CUS_URL\n",
    "\n",
    "(params['minElevation'], params['maxElevation']) = e.value\n",
    "params['output'] = 'PNG'\n",
    "params['peaks'] = peaks.value\n",
    "\n",
    "if c_sel.value == '--' and c2_sel.value != '--':\n",
    "    raise ValueError('Must specify primary to specify secondary')\n",
    "\n",
    "if g_sel.value == '--' and g2_sel.value != '--':\n",
    "    raise ValueError('Must specify primary to specify secondary')\n",
    "\n",
    "if c_sel.value != '--':\n",
    "    params['canopy_ds'] = c_sel.value\n",
    "    if c2_sel.value != '--':\n",
    "        params['canopy_ds'] = params['canopy_ds'] + ',' + c2_sel.value\n",
    "\n",
    "if g_sel.value != '--':\n",
    "    params['ground_ds'] = g_sel.value\n",
    "    if g2_sel.value != '--':\n",
    "        params['ground_ds'] = params['ground_ds'] + ',' + g2_sel.value\n",
    "\n",
    "params['elevPercentiles'] = percentile.value\n",
    "\n",
    "response = get(url, {p: params[p] for p in params if p not in ['bounds', 'kmzpng']})\n",
    "\n",
    "buf = BytesIO(response.content)\n",
    "\n",
    "img = Image.open(buf).convert('RGB')\n",
    "display(img)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d2d003c-a487-4747-81d0-bfcf88b32eca",
   "metadata": {},
   "source": [
    "## Tomogram Profile Plots - Fixed Examples"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca398588-8893-4ca4-a896-11d3a6448207",
   "metadata": {},
   "source": [
    "### Tomograms - Slicing along single lats/lons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b15334c-47a4-4fd5-8eb4-97dbe09e2e3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = dict(\n",
    "    ds=TOMOGRAM_DS_KC + '_NO_DEM',\n",
    "    output='PNG',\n",
    "    latitude=37.055,\n",
    "    minLon=-119.24,\n",
    "    maxLon=-119.2,\n",
    "    minElevation=-121,\n",
    "    maxElevation=121,\n",
    "    peaks=True\n",
    ")\n",
    "\n",
    "response = get(TOMOGRAM_LAT_URL, params)\n",
    "\n",
    "buf = BytesIO(response.content)\n",
    "\n",
    "img = Image.open(buf).convert('RGB')\n",
    "display(img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b41c622-5ba7-418c-a863-3e5a7b641001",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = dict(\n",
    "    ds=TOMOGRAM_DS_KC,\n",
    "    output='PNG',\n",
    "    longitude=-119.2103533733451,\n",
    "    minLat=37,\n",
    "    maxLat=37.25,\n",
    "    minElevation=1000,\n",
    "    maxElevation=2000,\n",
    "    peaks=False\n",
    ")\n",
    "\n",
    "response = get(TOMOGRAM_LON_URL, params)\n",
    "\n",
    "buf = BytesIO(response.content)\n",
    "\n",
    "img = Image.open(buf).convert('RGB')\n",
    "display(img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55f3c1b2-4366-496d-9444-c8235c1274ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = dict(\n",
    "    ds=TOMOGRAM_DS_Lope.format('VV'),\n",
    "    output='PNG',\n",
    "    longitude=11.6,\n",
    "    minLat=-0.075,\n",
    "    maxLat=-0.05,\n",
    "    minElevation=5,\n",
    "    maxElevation=350,\n",
    "    peaks=False\n",
    ")\n",
    "\n",
    "response = get(TOMOGRAM_LON_URL, params)\n",
    "\n",
    "buf = BytesIO(response.content)\n",
    "\n",
    "img = Image.open(buf).convert('RGB')\n",
    "display(img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efd9a26d-cc2d-4b12-aa6d-6fb583d765f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = dict(\n",
    "    ds=TOMOGRAM_DS_Rabi.format('HV'),\n",
    "    output='PNG',\n",
    "    latitude=-1.85,\n",
    "    minLon=9.8,\n",
    "    maxLon=10,\n",
    "    minElevation=-90,\n",
    "    maxElevation=1200,\n",
    "    peaks=False\n",
    ")\n",
    "\n",
    "response = get(TOMOGRAM_LAT_URL, params)\n",
    "\n",
    "buf = BytesIO(response.content)\n",
    "\n",
    "img = Image.open(buf).convert('RGB')\n",
    "display(img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "834699ca-98dd-4363-8c7f-a71cb5f54541",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = dict(\n",
    "    ds=TOMOGRAM_DS_Rabi.format('HV') + '_NO_DEM',\n",
    "    output='PNG',\n",
    "    longitude=9.98,\n",
    "    minLat=-1.972,\n",
    "    maxLat=-1.935,\n",
    "    minElevation=-40,\n",
    "    maxElevation=120,\n",
    "    peaks=False,\n",
    "    canopy_ds='rabi_L_KauriCNN_canopy_height_5x5_200epochs',\n",
    "    ground_ds='rabi_L_KauriCNN_ground_elevation_bias_5x5_60epochs',\n",
    "    elevPercentiles=True\n",
    ")\n",
    "\n",
    "response = get(TOMOGRAM_LON_URL, params)\n",
    "\n",
    "buf = BytesIO(response.content)\n",
    "\n",
    "img = Image.open(buf).convert('RGB')\n",
    "display(img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77919217-ad6f-403e-be09-6802ded3236b",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = dict(\n",
    "    ds=TOMOGRAM_DS_Rabi.format('HV') + '_NO_DEM',\n",
    "    output='PNG',\n",
    "    longitude=9.98,\n",
    "    minLat=-1.972,\n",
    "    maxLat=-1.935,\n",
    "    minElevation=-40,\n",
    "    maxElevation=120,\n",
    "    peaks=False,\n",
    "    canopy_ds='rabi_L_KauriCNN_canopy_height_5x5_200epochs',\n",
    "    ground_ds='rabi_L_KauriCNN_ground_elevation_bias_5x5_60epochs',\n",
    "    elevPercentiles=True,\n",
    "    cmap='forest'\n",
    ")\n",
    "\n",
    "response = get(TOMOGRAM_LON_URL, params)\n",
    "\n",
    "buf = BytesIO(response.content)\n",
    "\n",
    "img = Image.open(buf).convert('RGB')\n",
    "display(img)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72e31def-dc10-4260-9e9d-cc94c4a43488",
   "metadata": {},
   "source": [
    "### Tomograms - Slicing along range\n",
    "<!-- is there a better way to title this section??? /-->\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7c9b789-1cc8-4bbf-a60d-ee9f3c903bd8",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = dict(\n",
    "    ds=TOMOGRAM_DS_KC,\n",
    "    output='GIF',\n",
    "    longitude='-119.2103533733451:-119.2082533733451',\n",
    "    minLat=37,\n",
    "    maxLat=37.25,\n",
    "    minElevation=1000,\n",
    "    maxElevation=2000,\n",
    "    peaks=True,\n",
    ")\n",
    "\n",
    "response = get(TOMOGRAM_LON_URL, params)\n",
    "\n",
    "img = IpImage(data=response.content, format='gif')\n",
    "display(img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eff9ba4e-31f3-45ff-ab28-5e86996130ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = dict(\n",
    "    ds=TOMOGRAM_DS_Rabi.format('HV'),\n",
    "    output='GIF',\n",
    "    longitude='9.97:9.99',\n",
    "    minLat=-1.97,\n",
    "    maxLat=-1.93,\n",
    "    minElevation=10,\n",
    "    maxElevation=150,\n",
    "    peaks=False,\n",
    "    canopy_ds='rabi_L_KauriCNN_canopy_elevation_5x5_200epochs',\n",
    "    ground_ds='rabi_L_KauriCNN_ground_elevation_dem_5x5_60epochs'\n",
    ")\n",
    "\n",
    "response = get(TOMOGRAM_LON_URL, params)\n",
    "\n",
    "img = IpImage(data=response.content, format='gif')\n",
    "display(img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a94971e4-d022-487c-9210-caf29cb5dd4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = dict(\n",
    "    ds=TOMOGRAM_DS_Rabi.format('HV'),\n",
    "    output='GIF',\n",
    "    longitude='9.97:9.99',\n",
    "    minLat=-1.97,\n",
    "    maxLat=-1.93,\n",
    "    minElevation=10,\n",
    "    maxElevation=150,\n",
    "    peaks=False,\n",
    "    canopy_ds='rabi_L_KauriCNN_canopy_elevation_5x5_200epochs',\n",
    "    ground_ds='rabi_L_KauriCNN_ground_elevation_dem_5x5_60epochs',\n",
    "    elevPercentiles=True\n",
    ")\n",
    "\n",
    "response = get(TOMOGRAM_LON_URL, params)\n",
    "\n",
    "img = IpImage(data=response.content, format='gif')\n",
    "display(img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f0999e9-5672-41c4-bf0e-09b9e40cef41",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = dict(\n",
    "    ds=TOMOGRAM_DS_Rabi.format('HH'),\n",
    "    output='GIF',\n",
    "    latitude='-1.86:-1.84',\n",
    "    minLon=9.8,\n",
    "    maxLon=9.935,\n",
    "    minElevation=-90,\n",
    "    maxElevation=1200,\n",
    "    peaks=False,\n",
    ")\n",
    "\n",
    "response = get(TOMOGRAM_LAT_URL, params)\n",
    "\n",
    "img = IpImage(data=response.content, format='gif')\n",
    "display(img)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c19fb1cb-51bf-46b8-890b-5ee23f52b3f4",
   "metadata": {},
   "source": [
    "### Tomograms - 3D Plot Examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c53a05aa-3365-4094-a478-10f7c1e6bfd8",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = dict(\n",
    "    ds=TOMOGRAM_DS_Rabi.format('VV'),\n",
    "    output='PNG',\n",
    "    b='9.85,-1.87,9.86,-1.86',\n",
    "    minElevation=-95,\n",
    "    maxElevation=240,\n",
    "    viewAngle='225,15'\n",
    ")\n",
    "\n",
    "response = get(TOMOGRAM_3D_URL, params)\n",
    "\n",
    "buf = BytesIO(response.content)\n",
    "\n",
    "img = Image.open(buf).convert('RGB')\n",
    "display(img)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "305f6286-f2d3-4ede-8697-f0e9d3f226bf",
   "metadata": {},
   "source": [
    "#### 3D Masking & Rendering options"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5acb9aa9-c962-4062-8e6d-f24b22ab412a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Raw tomogram\n",
    "\n",
    "params = dict(\n",
    "    ds=TOMOGRAM_DS_Rabi.format('VV') + '_NO_DEM',\n",
    "    output='PNG',\n",
    "    b='9.95,-1.943,10,-1.893',\n",
    "    minElevation=-102,\n",
    "    maxElevation=56,\n",
    "    viewAngle='225,15'\n",
    ")\n",
    "\n",
    "response = get(TOMOGRAM_3D_URL, params)\n",
    "\n",
    "buf = BytesIO(response.content)\n",
    "\n",
    "img = Image.open(buf).convert('RGB')\n",
    "display(img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99b8be6d-f53e-437e-b183-a55d3fadeb94",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Masked tomogram\n",
    "\n",
    "params = dict(\n",
    "    ds=TOMOGRAM_DS_Rabi.format('VV') + '_NO_DEM',\n",
    "    output='PNG',\n",
    "    b='9.95,-1.943,10,-1.893',\n",
    "    minElevation=-40,\n",
    "    maxElevation=56,\n",
    "    viewAngle='225,30',\n",
    "    canopy_ds='rabi_L_KauriCNN_canopy_height_5x5_200epochs',\n",
    "    ground_ds='rabi_L_KauriCNN_ground_elevation_bias_5x5_60epochs',\n",
    ")\n",
    "\n",
    "response = get(TOMOGRAM_3D_URL, params)\n",
    "\n",
    "buf = BytesIO(response.content)\n",
    "\n",
    "img = Image.open(buf).convert('RGB')\n",
    "display(img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96cc5408-76af-4a27-bd81-a2197dbbb11a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Percentile rendering\n",
    "\n",
    "params = dict(\n",
    "    ds=TOMOGRAM_DS_Rabi.format('VV') + '_NO_DEM',\n",
    "    output='PNG',\n",
    "    b='9.95,-1.943,10,-1.893',\n",
    "    minElevation=-40,\n",
    "    maxElevation=56,\n",
    "    viewAngle='225,30',\n",
    "    canopy_ds='rabi_L_KauriCNN_canopy_height_5x5_200epochs',\n",
    "    ground_ds='rabi_L_KauriCNN_ground_elevation_bias_5x5_60epochs',\n",
    "    elevPercentiles=True,\n",
    ")\n",
    "\n",
    "response = get(TOMOGRAM_3D_URL, params)\n",
    "\n",
    "buf = BytesIO(response.content)\n",
    "\n",
    "img = Image.open(buf).convert('RGB')\n",
    "display(img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e6cec490-f451-443b-9340-f0a20fedc7f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# RH98 rendering w/ custom color map\n",
    "\n",
    "params = dict(\n",
    "    ds=TOMOGRAM_DS_Rabi.format('VV') + '_NO_DEM',\n",
    "    output='PNG',\n",
    "    b='9.95,-1.943,10,-1.893',\n",
    "    minElevation=-40,\n",
    "    maxElevation=56,\n",
    "    viewAngle='225,30',\n",
    "    canopy_ds='rabi_L_KauriCNN_canopy_height_5x5_200epochs',\n",
    "    ground_ds='rabi_L_KauriCNN_ground_elevation_bias_5x5_60epochs',\n",
    "    elevPercentiles=True,\n",
    "    renderRH=98,\n",
    "    cmap='forest',\n",
    ")\n",
    "\n",
    "response = get(TOMOGRAM_3D_URL, params)\n",
    "\n",
    "buf = BytesIO(response.content)\n",
    "\n",
    "img = Image.open(buf).convert('RGB')\n",
    "display(img)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4bfeba35-5d02-4d4a-addf-2c44f8e0a742",
   "metadata": {},
   "outputs": [],
   "source": [
    "# RH98 rendering w/ custom color map\n",
    "\n",
    "params = dict(\n",
    "    ds=TOMOGRAM_DS_Rabi.format('VV') + '_NO_DEM',\n",
    "    output='PNG',\n",
    "    b='9.97,-1.943,9.99,-1.923',\n",
    "    minElevation=-40,\n",
    "    maxElevation=56,\n",
    "    viewAngle='225,30',\n",
    "    canopy_ds='rabi_L_KauriCNN_canopy_height_5x5_200epochs',\n",
    "    ground_ds='rabi_L_KauriCNN_ground_elevation_bias_5x5_60epochs',\n",
    "    elevPercentiles=True,\n",
    "    renderRH=98,\n",
    "    cmap='forest',\n",
    ")\n",
    "\n",
    "response = get(TOMOGRAM_3D_URL, params)\n",
    "\n",
    "buf = BytesIO(response.content)\n",
    "\n",
    "img = Image.open(buf).convert('RGB')\n",
    "display(img)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef633b44-61da-4c33-9e75-ed41c55ee39b",
   "metadata": {},
   "source": [
    "### Tomograms - Other Output Formats\n",
    "\n",
    "Tomograms subsets can be returned in a simple CSV format to be used with VR tools to enable users to interactively manipulate and explore the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e684270-9994-4eab-9d10-cde7a0c11cb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = dict(\n",
    "    ds=TOMOGRAM_DS_Rabi.format('VV'),\n",
    "    output='CSV',\n",
    "    b='9.95,-1.97,10,-1.92',\n",
    "    minElevation=-95,\n",
    "    maxElevation=240\n",
    ")\n",
    "\n",
    "response = get(TOMOGRAM_3D_URL, params)\n",
    "\n",
    "buf = StringIO(response.text)\n",
    "\n",
    "df = pd.read_csv(buf, names=['x', 'y', 'z', 'value'])\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d4b38a8-e914-46ae-ae50-56d34ea9f8c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "df.to_csv('tomo.ply', header=False, index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aff0c8da-3d5c-44bc-b45b-c02ab5329921",
   "metadata": {},
   "source": [
    "## LIDAR\n",
    "\n",
    "The following cells will showcase SDAP's LIDAR endpoint (Still a WIP)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "52cc3641-3b53-4223-8def-ff21af06d50c",
   "metadata": {},
   "source": [
    "### LIDAR - PNG Output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "145edd7e-7065-4012-88ad-2cd796590617",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = dict(\n",
    "    ds=LIDAR_DS,\n",
    "    b='-148.8,66.025,-148.5,66.18',\n",
    "    startTime='2016-06-28T00:00:00Z',\n",
    "    endTime='2020-06-30T00:00:00Z',\n",
    "    output='PNG',\n",
    "    mapToGrid=True\n",
    ")\n",
    "\n",
    "response = get(LIDAR_URL, params)\n",
    "response.raise_for_status()\n",
    "buf = BytesIO(response.content)\n",
    "\n",
    "img = Image.open(buf).convert('RGB')\n",
    "display(img)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40c93d3d-9c11-47ba-a00d-ee061007312b",
   "metadata": {},
   "source": [
    "### LIDAR - NetCDF Output & Profile Slicing\n",
    "\n",
    "Profile slicing can be defined by longitude, latitude, or an arbitrary WKT LineString.\n",
    "\n",
    "It is supported by all output formats except CSV, but showcased here for the NetCDF output format as the plotting in the PNG format is still an awkward work in progress."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd680957",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = dict(\n",
    "    ds=LIDAR_DS,\n",
    "    b='-122.3,46,-122,46.25',\n",
    "    startTime='2016-06-28T00:00:00Z',\n",
    "    endTime='2020-06-30T00:00:00Z',\n",
    "    output='NetCDF',\n",
    "    filename='lidar.nc',\n",
    "    mapToGrid=True,\n",
    "    sliceSamples=10000,\n",
    "    sliceWKT='LINESTRING (-122.12 46.07, -122.09095182261792 46.0220916078604)',\n",
    "    latSlice=46.09\n",
    ")\n",
    "\n",
    "response = get(LIDAR_URL, params)\n",
    "response.raise_for_status()\n",
    "buf = BytesIO(response.content)\n",
    "\n",
    "ds = xr.open_dataset(buf).isel(time=0)\n",
    "ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6cac979-9002-4488-99e3-0380c4e30ab5",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.ground_height.plot(cmap='viridis'); plt.show()\n",
    "ds.mean_veg_height.plot(cmap='viridis', vmin=0); plt.show()\n",
    "ds.canopy_height.plot(cmap='viridis', vmin=0, vmax=60); plt.show()\n",
    "ds.canopy_coverage.plot(cmap='viridis', vmin=0, vmax=100); plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0dd9762-99cb-49a2-b8a5-8afbcf52ab15",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_pts = ds.lon.to_numpy()\n",
    "rh50_pts = ds['lat_slice_mean_veg_height'].to_numpy()\n",
    "rh98_pts = ds['lat_slice_canopy_height'].to_numpy()\n",
    "zg_pts = ds['lat_slice_ground_height'].to_numpy()\n",
    "cc_pts = ds['lat_slice_canopy_coverage'].to_numpy()\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(2,1, figsize=(10, 8))\n",
    "\n",
    "ax1.plot(\n",
    "    x_pts, zg_pts,\n",
    "    x_pts, zg_pts + rh50_pts,\n",
    "    x_pts, zg_pts + rh98_pts,\n",
    ")\n",
    "ax1.set_title('Longitudunal profile heights (m)')\n",
    "ax1.legend([\n",
    "    'Ground Height',\n",
    "    'Mean Vegetation Height',\n",
    "    'Canopy Height'\n",
    "])\n",
    "ax1.set_ylabel('Height [m]')\n",
    "\n",
    "ax2.plot(x_pts, cc_pts)\n",
    "ax2.set_title('Longitudunal profile canopy coverage (%)')\n",
    "ax2.legend(['Canopy Coverage'])\n",
    "ax2.set_ylim([0, 100])\n",
    "ax2.set_xlabel('Longitude [deg]')\n",
    "ax2.set_ylabel('Height [m]')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "946daa00-31d5-45c8-b1e4-fad7bbd12601",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_pts = ds.distance.to_numpy()\n",
    "rh50_pts = ds['slice_mean_veg_height'].to_numpy()\n",
    "rh98_pts = ds['slice_canopy_height'].to_numpy()\n",
    "zg_pts = ds['slice_ground_height'].to_numpy()\n",
    "cc_pts = ds['slice_canopy_coverage'].to_numpy()\n",
    "\n",
    "fig, (ax1, ax2) = plt.subplots(2,1, figsize=(10, 8))\n",
    "\n",
    "ax1.plot(\n",
    "    x_pts, zg_pts,\n",
    "    x_pts, zg_pts + rh50_pts,\n",
    "    x_pts, zg_pts + rh98_pts,\n",
    ")\n",
    "ax1.set_title('Custom slice heights (m)')\n",
    "ax1.legend([\n",
    "    'Ground Height',\n",
    "    'Mean Vegetation Height',\n",
    "    'Canopy Height'\n",
    "])\n",
    "ax1.set_ylabel('Height [m]')\n",
    "\n",
    "ax2.plot(x_pts, cc_pts)\n",
    "ax2.set_title('Custom slice canopy coverage (%)')\n",
    "ax2.legend(['Canopy Coverage'])\n",
    "ax2.set_ylim([0, 100])\n",
    "ax2.set_ylabel('Height [m]')\n",
    "ax2.set_xlabel('Distance along line [m]')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eed0c64b-684b-4fe5-afdc-09e2694b4f0e",
   "metadata": {},
   "source": [
    "### LIDAR - CSV & ZIP Outputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "718e6232-ebcf-4da7-80af-1baca10f462c",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = dict(\n",
    "    ds=LIDAR_DS,\n",
    "    b='-122.3,46,-122,46.25',\n",
    "    startTime='2016-06-28T00:00:00Z',\n",
    "    endTime='2020-06-30T00:00:00Z',\n",
    "    output='CSV',\n",
    "    filename='lidar.csv',\n",
    "    mapToGrid=True,\n",
    ")\n",
    "\n",
    "response = get(LIDAR_URL, params)\n",
    "response.raise_for_status()\n",
    "buf = StringIO(response.text)\n",
    "\n",
    "df = pd.read_csv(buf)\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "beedf44a-586d-4b18-9dbf-c7a1609a204e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# .ply formatting\n",
    "\n",
    "df[['longitude', 'latitude', 'canopy_height', 'canopy_coverage']].to_csv('lidar.ply', header=False, index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9b9b248-e109-420f-96a0-3bc6173ccf35",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = dict(\n",
    "    ds=LIDAR_DS,\n",
    "    b='-122.3,46,-122,46.25',\n",
    "    startTime='2016-06-28T00:00:00Z',\n",
    "    endTime='2020-06-30T00:00:00Z',\n",
    "    output='ZIP',\n",
    "    filename='lidar.zip',\n",
    "    mapToGrid=True,\n",
    "    sliceSamples=10000,\n",
    "    sliceWKT='LINESTRING (-122.12 46.07, -122.09095182261792 46.0220916078604)',\n",
    "    latSlice='46.09,-122.15,-122.11'\n",
    ")\n",
    "\n",
    "response = get(LIDAR_URL, params)\n",
    "response.raise_for_status()\n",
    "buf = BytesIO(response.content)\n",
    "\n",
    "zipfile = ZipFile(buf)\n",
    "zipfile.namelist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2227e86c-2ce0-41dd-95d0-d81cdadae842",
   "metadata": {},
   "outputs": [],
   "source": [
    "for csv in zipfile.namelist():\n",
    "    print(f'{csv}:')\n",
    "    csv_buf = StringIO(zipfile.read(csv).decode('utf-8'))\n",
    "    csv_data = pd.read_csv(csv_buf)\n",
    "    display(csv_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "968845a0-fe58-4ee0-bd58-e0e31c8b6e5d",
   "metadata": {},
   "source": [
    "## Kauri output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee6011c4-2368-437f-adf7-d75f654fa85b",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_epochs_map = {  # Map number of training epochs used for the model to the timestamp of the data product it created\n",
    "    20: '2024-04-16T11:12:00Z',\n",
    "    60: '2024-05-07T17:54:36Z',\n",
    "    200: '2024-05-07T22:54:05Z' if not LOCAL else '2024-05-07T17:54:36Z'\n",
    "}\n",
    "\n",
    "bbox_map = dict(\n",
    "    full='9.70,-2.05,10.11,-1.65',\n",
    "    closeup='9.95,-2.0,10.05,-1.9'\n",
    ")\n",
    "\n",
    "output_map = {}\n",
    "\n",
    "def kauri(b, n_epochs, ax):\n",
    "    if n_epochs not in n_epochs_map:\n",
    "        raise ValueError(f'No associated timestamp for {n_epochs=}')\n",
    "\n",
    "    timestamp = n_epochs_map[n_epochs]\n",
    "    \n",
    "    params = {\n",
    "        'dataset' : KAURI_SAMPLE_DS,\n",
    "        'b' : bbox_map[b],\n",
    "        'startTime' : timestamp, \n",
    "        'endTime' : timestamp,\n",
    "        'output' : 'ZIP'\n",
    "    }\n",
    "    \n",
    "    url = f'{NEXUS_URL}/cdmssubset'\n",
    "    \n",
    "    t = datetime.now()\n",
    "        \n",
    "    try:\n",
    "        print('Waiting for response from SDAP...', end='')\n",
    "        response = requests.get(url, params=params, verify=False)\n",
    "        td = datetime.now() - t\n",
    "        response.raise_for_status()\n",
    "    \n",
    "        print(f' Done in {td}')\n",
    "    except:\n",
    "        td = datetime.now() - t\n",
    "        print(f' Failed in {td}')\n",
    "        print(response.text)\n",
    "        raise\n",
    "        \n",
    "    response = BytesIO(response.content)\n",
    "    \n",
    "    print('Unpacking results')\n",
    "    \n",
    "    with ZipFile(response) as zipfile:\n",
    "        namelist = zipfile.namelist()\n",
    "        csv_buf = StringIO(zipfile.read(namelist[0]).decode('utf-8'))\n",
    "        csv_data = pd.read_csv(csv_buf)\n",
    "\n",
    "        output_map[(b, n_epochs)] = csv_data\n",
    "        \n",
    "    print('Getting coordinate arrays')\n",
    "    lats = np.unique(csv_data['latitude'].values)\n",
    "    lons = np.unique(csv_data['longitude'].values)\n",
    "    times = np.unique(csv_data['time'].values)\n",
    "    \n",
    "    data_col = [c for c in csv_data.columns if c not in ['time', 'latitude', 'longitude']][0]\n",
    "    \n",
    "    print('Gridding the result')\n",
    "    \n",
    "    vals_3d = np.empty((len(times), len(lats), len(lons)), dtype='float32')\n",
    "    data_dict = {}\n",
    "    \n",
    "    for t in csv_data.itertuples(index=False):\n",
    "        key = (t.time, t.latitude, t.longitude)\n",
    "        data_dict[key] = getattr(t, data_col)\n",
    "    \n",
    "    for i, t in enumerate(times):\n",
    "        for j, lat in enumerate(lats):\n",
    "            for k, lon in enumerate(lons):\n",
    "                vals_3d[i, j, k] = data_dict.get((t, lat, lon), np.nan)\n",
    "            \n",
    "\n",
    "    v_min = np.nanmin(vals_3d)\n",
    "    v_max = np.nanmax(vals_3d)\n",
    "    \n",
    "    pcm = ax.pcolormesh(*np.meshgrid(lons, lats), vals_3d[0])\n",
    "    ax.set_xlabel('Longitude')\n",
    "    ax.set_ylabel('Latitude')\n",
    "    ax.set_title(f'{n_epochs} epochs: {b}')\n",
    "    return [v_min, v_max]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e45f2186-d5de-4050-b32b-432e7702b0f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 2, figsize=(14,12))\n",
    "\n",
    "minmax = []\n",
    "\n",
    "minmax.append(kauri('full', 20, axes[0][0]))\n",
    "minmax.append(kauri('closeup', 20, axes[0][1]))\n",
    "minmax.append(kauri('full', 200, axes[1][0]))\n",
    "minmax.append(kauri('closeup', 200, axes[1][1]))\n",
    "\n",
    "minmax = np.array(minmax)\n",
    "\n",
    "cb = fig.colorbar(\n",
    "    ScalarMappable(norm=Normalize(vmin=min(minmax[:, 0]), vmax=max(minmax[:, 1])), cmap='viridis'), \n",
    "    label='Canopy height (98th percentile) [m]', \n",
    "    ax=axes.ravel().tolist()\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b89a562-2e22-4a87-bb2a-e897563ff8c2",
   "metadata": {},
   "outputs": [],
   "source": [
    "for o in output_map:\n",
    "    output = output_map[o].drop('time', axis=1)\n",
    "    col = output.columns[-1]\n",
    "    output['val'] = output[col]\n",
    "    output.to_csv(f'kauri_{o[0]}_{o[1]}.ply', header=False, index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb9a7017-e136-49c2-a1d5-e49b82dd0f85",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:stv-fis] *",
   "language": "python",
   "name": "conda-env-stv-fis-py"
  },
  "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}