is2_quicklook_depths.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6b8bdc6f-133f-4abd-b0e2-51d75c7ac4d5",
   "metadata": {
    "user_expressions": []
   },
   "source": [
    "# Estimating Snow Depth from ICESat-2 ATL08QL data\n",
    "\n",
    "This script is made to derive preliminary snow depth estimates using the ATL08 Quick Look product. Because the quick looks are (a) subject to geolocation errors and (b) automatically deleted after the final data release, it is recommended to use `is2_snowex_ak.ipynb` if it has been ~45 days since the RGT of interest.\n",
    "\n",
    "icepyx normally does not support downloading of ATL08QL data, so a slightly altered version is required. The user may find a branch with ATL08QL functionality at: https://github.com/zachghiaccio/icepyx/tree/quick-looks. A pull request for this feature will be introduced to the icepyx development team in the near future.\n",
    "\n",
    "When ATL08QL functionality is introduced to icepyx, it is planned to integrate this notebook with `is2_snowex_ak.ipynb`.\n",
    "\n",
    "As with the above notebook, snow-on/-off DEMs from the UAF airborne lidar are needed to run this script. Interested users should contact the original provider of the data (Chris Larsen, cflarsen@uaf.edu)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7b5a787-c646-4d94-9c15-967942c97348",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import h5py\n",
    "import matplotlib.pyplot as plt\n",
    "import rioxarray as rio\n",
    "import os\n",
    "from osgeo import gdal\n",
    "import icepyx as ipx\n",
    "import pandas as pd\n",
    "import pickle\n",
    "from pyproj import Proj, transform, Transformer, CRS\n",
    "import sys\n",
    "\n",
    "sys.path.append('C:/Users/zfair/OneDrive - NASA/Documents/Python/icesat2-snow/')\n",
    "import lidar_processing as lp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "473c2763-588f-4e81-9e8f-7528383c6548",
   "metadata": {
    "user_expressions": []
   },
   "source": [
    "To keep things simple, we are not using cloud services to obtain the ATL08QL data - the data will be downloaded to our workspace. Once this notebook is integrated, the hope is to allow cloud access to quick look data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2e3eaa1-63b6-42f2-9667-a4d832ee0438",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify the ICESat-2 product\n",
    "short_name = 'ATL08QL'\n",
    "\n",
    "# Define the spatial extent using a pre-generated bounding box\n",
    "with open('C:/Users/zfair/OneDrive - NASA/Documents/Python/icesat2-snow/Data/shapefiles/snowex_sites_for_icepyx.pkl', 'rb') as f:\n",
    "    coordinates = pickle.load(f)\n",
    "    spatial_extent = coordinates['alaska']\n",
    "    \n",
    "# Date range of requested data\n",
    "date_range = ['2023-03-01', '2023-03-31']\n",
    "\n",
    "# ICESat-2 track (RGT, optional)\n",
    "rgt = '1356'\n",
    "    \n",
    "# Generate the query object\n",
    "try:\n",
    "    # With RGT specification\n",
    "    region = ipx.Query(short_name, spatial_extent, date_range, tracks=rgt)\n",
    "except:\n",
    "    # General query\n",
    "    region = ipx.Query(short_name, spatial_extent, date_range)\n",
    "    \n",
    "print(region)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "167817e5-261a-43ae-973c-07f18ffce758",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Enter login information for Earthdata\n",
    "uid = 'zhfair'\n",
    "email = 'zhfair@umich.edu'\n",
    "region.earthdata_login(uid, email)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0f5b802",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check the ids for the requested data\n",
    "region.avail_granules(ids=True)\n",
    "\n",
    "date_str = region.avail_granules(ids=True)[0][0][8:16]\n",
    "print(date_str)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d73f7dc-c5fb-42e1-bbfa-b645f030550d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Start the data order\n",
    "region.order_granules()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb5544ef-c125-468f-a163-eed6f22c7d0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download the granules into the specified folder\n",
    "path = 'C:/Users/zfair/OneDrive - NASA/Documents/Python/icesat2-snow/Data/is2/icepyx/'\n",
    "region.download_granules(path)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19e1576c-25a5-483a-8275-9f19d0a1ad38",
   "metadata": {
    "user_expressions": []
   },
   "source": [
    "Now that our data is downloaded, let's load the file into a DataFrame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c513327-25da-4f7a-9545-f0b882c81525",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find file in saved directory\n",
    "for fname in os.listdir(path):\n",
    "    if rgt and date_str in fname:\n",
    "        is2_file = path + fname\n",
    "                    \n",
    "print(is2_file)\n",
    "                    \n",
    "# Load the ATL08QL data into a DataFrame\n",
    "with h5py.File(is2_file, 'r') as f:\n",
    "    # Identify the strong beams\n",
    "    sc_orient = f['orbit_info/sc_orient'][0]\n",
    "    strong_beams, strong_ids = lp.strong_beam_finder(sc_orient)\n",
    "    \n",
    "    # Concatenate strong beam data into DataFrame\n",
    "    atl08 = lp.beam_cycle_concat([is2_file], 'ATL08')\n",
    "    \n",
    "atl08.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02c76774",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add easting/northing coordinates (EPSG:32606)\n",
    "transformer = Transformer.from_crs('EPSG:4326', 'EPSG:32606', always_xy=True)\n",
    "\n",
    "atl08['x'], atl08['y'] = transformer.transform(atl08.lon, atl08.lat)\n",
    "\n",
    "atl08.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1535575c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove messy ICESat-2 data\n",
    "upper = atl08.height.mean() + 3*atl08.height.std()\n",
    "atl08 = atl08.loc[atl08.height<upper]\n",
    "\n",
    "# Add vertical height correction needed for quick-look products\n",
    "atl08['height'] += 2.7"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "407d0622",
   "metadata": {},
   "source": [
    "Our ICESat-2 data is now ready! Let's load up a snow-off DEM from the UAF airborne lidar."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2777b3fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load UAF lidar data into a raster. Will be updated to be more flexible at a later date.\n",
    "f_snow_off = 'C:/Users/zfair/OneDrive - NASA/Documents/Python/icesat2-snow/Data/UAF/farmersloop_2022may28_dtm.tif'\n",
    "f_snow_on = 'C:/Users/zfair/OneDrive - NASA/Documents/Python/icesat2-snow/Data/UAF/farmersloop_2022oct24_snowdepth.tif'\n",
    "\n",
    "lidar_snow_off = rio.open_rasterio(f_snow_off)\n",
    "lidar_snow_on = rio.open_rasterio(f_snow_on)\n",
    "\n",
    "lidar_snow_off"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1758d355",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply vertical datum correction. To be updated to be more flexible.\n",
    "lidar_snow_off -= 9.75\n",
    "\n",
    "# Coregister with ICESat-2 using a bivariate spline\n",
    "strong_ids = np.unique(atl08['gt'])\n",
    "atl08_uaf = lp.coregister_is2(lidar_snow_off, lidar_snow_on, atl08, strong_ids)\n",
    "\n",
    "atl08_uaf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50a47267",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove filler data\n",
    "atl08_uaf['lidar_height'][atl08_uaf['lidar_height']<0] = np.nan\n",
    "atl08_uaf['lidar_snow_depth'][atl08_uaf['lidar_snow_depth']<0] = np.nan\n",
    "atl08_uaf['residual'][atl08_uaf['residual']>100] = np.nan\n",
    "\n",
    "atl08_uaf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8eb863e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "atl08_uaf.plot.scatter(x='y', y='residual')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4dc17e9d",
   "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.10.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}