ICESat-2_DEM_comparison_Colombia_working.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparing ICESat-2 Altimetry Elevations with DEM\n",
    "### Example Notebook\n",
    "This notebook compares elevations from ICESat-2 to those from a DEM.\n",
    "\n",
    "#### Credits\n",
    "* notebook by: [Jessica Scheick](https://github.com/JessicaS11) and [Shashank Bhushan](https://github.com/ShashankBice)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Setup\n",
    "\n",
    "If you are not running this notebook on the [Icepyx Gallery Image](https://hub.docker.com/r/icepyximg/icepyx-gallery-image), please use the following to set up your environment.\n",
    "\n",
    "##### The Notebook was run on ICESat2 Hackweek 2019 pangeo image\n",
    "\n",
    "##### For full functionality,\n",
    "\n",
    "- Please install [icepyx](https://github.com/icesat2py/icepyx), [topolib](https://github.com/ICESAT-2HackWeek/topohack), [contextily](https://github.com/darribas/contextily) using `git clone xxxxx`, `pip install -e .` workflow (see below; **you must restart your kernel after installing the packages**)\n",
    "- Download [NASA ASP](https://github.com/NeoGeographyToolkit/StereoPipeline) tar ball and unzip, we execute the commands from the notebook, using the path to the untared bin folder for the given commands."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "%%bash\n",
    "cd ~\n",
    "# git clone https://github.com/icesat2py/icepyx.git\n",
    "# git clone https://github.com/ICESAT-2HackWeek/topohack.git\n",
    "# git clone https://github.com/darribas/contextily.git\n",
    "\n",
    "cd contextily\n",
    "pip install -e .\n",
    "cd ../topohack\n",
    "pip install -e .\n",
    "cd ../icepyx\n",
    "pip install -e .\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "%cd ~\n",
    "#needs to be wherever icepyx, contextily, and topolib are installed in the previous step (ideally $HOME)\n",
    "# %pwd\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you are using the Icepyx Gallery Image, continue on"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### ICESat-2 product being explored : [ATL08](https://nsidc.org/data/atl08)\n",
    "- Along track heights for canopy (land and vegitation) and  terrain\n",
    "- Terrain heights provided are aggregated over every 100 m along track interval, output contains \"h_te_best_fit: height from best fit algorithm for all photons in the range\", median height and others. Here we use h_te_best_fit.\n",
    "- See this preliminary introduction and quality assessment [paper](https://www.mdpi.com/2072-4292/11/14/1721) for more detail"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import packages, including icepyx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import icepyx as ipx\n",
    "import os\n",
    "import shutil\n",
    "import h5py\n",
    "import xarray as xr\n",
    "# depedencies\n",
    "import getpass\n",
    "#from topolib.subsetDat import subsetBBox;\n",
    "from topolib import icesat2_data\n",
    "import glob\n",
    "import rasterio\n",
    "from topolib import gda_lib\n",
    "from topolib import dwnldArctic\n",
    "import numpy as np\n",
    "import geopandas as gpd\n",
    "from multiprocessing import Pool\n",
    "import contextily as ctx\n",
    "import pandas as pd\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%cd ~/icepyx/doc/examples/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preprocess #1\n",
    "- Download using icepyx"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Create an ICESat-2 data object with the desired search parameters\n",
    "- See the ICESat-2 DAAC Data Access notebook for more details on downloading data from the NSIDC"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "region_a = ipx.Query('ATL08', [-73.9, 10.7, -73.4, 11.1], ['2018-12-01','2019-09-01'], \\\n",
    "                          start_time='00:00:00', end_time='23:59:59')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finding and downloading data\n",
    "In order to download any data from NSIDC, we must first authenticate ourselves using a valid Earthdata login (available for free on their website). This will create a valid token to interface with the DAAC as well as start an active logged-in session to enable data download. The token is attached to the data object and stored, but the session must be passed to the download function. Then we can order the granules."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Log in to Earthdata"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "earthdata_uid = 'Jessica.scheick'\n",
    "email = 'jessica.scheick@maine.edu'\n",
    "region_a.earthdata_login(earthdata_uid, email)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#search for available granules\n",
    "region_a.avail_granules()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "region_a.granules.avail"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Place the order"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "region_a.order_granules(subset=False)\n",
    "#region_a.order_granules(verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#view a short list of order IDs\n",
    "region_a.granules.orderIDs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Download the order\n",
    "Finally, we can download our order to a specified directory (which needs to have a full path but doesn't have to point to an existing directory) and the download status will be printed as the program runs. Additional information is again available by using the optional boolean keyword 'verbose'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wd=%pwd\n",
    "path = wd + '/download'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "region_a.download_granules(path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Clean up the download folder by removing individual order folders:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Clean up Outputs folder by removing individual granule folders \n",
    "for root, dirs, files in os.walk(path, topdown=False):\n",
    "    for file in files:\n",
    "        try:\n",
    "            shutil.move(os.path.join(root, file), path)\n",
    "        except OSError:\n",
    "            pass\n",
    "        \n",
    "for root, dirs, files in os.walk(path):\n",
    "    for name in dirs:\n",
    "        os.rmdir(os.path.join(root, name))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preprocess #2\n",
    "- Convert data into geopandas dataframe, which allows for doing basing geospatial opertaions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# glob to list of files (run block of code creating wd and path variables if starting processing here)\n",
    "ATL08_list = sorted(glob.glob(path+'/*.h5'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Examine content of 1 ATLO8 hdf file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = ATL08_list[5]\n",
    "with h5py.File(filename, 'r') as f:\n",
    "    # List all groups\n",
    "    pairs=[1, 2, 3]\n",
    "    beams=['l','r']\n",
    "    print(\"Keys: %s\" % f.keys())\n",
    "    a_group_key = list(f.keys())[0]\n",
    "#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ATL08_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# dict containing data entries to retrive\n",
    "dataset_dict = {'land_segments':['delta_time','longitude','latitude','atl06_quality_summary','quality','terrain_flg'], 'land_segments/terrain':['h_te_best_fit']}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#gda_lib.ATL08_to_dict(ATL08_list[0],dataset_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## the data can be converted to geopandas dataframe, see ATL08_2_gdf function in topolib gda_lib\n",
    "temp_gdf = gda_lib.ATL08_2_gdf(ATL08_list[0],dataset_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temp_gdf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temp_gdf.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(temp_gdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colombia_crs = {'init':'epsg:32618'}\n",
    "plot_web = {'init':'epsg:3857'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "temp_gdf.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convert the list of hdf5 files into more familiar Pandas Dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_list = [(gda_lib.ATL08_2_gdf(x,dataset_dict)) for x in ATL08_list]\n",
    "gdf_colombia = gda_lib.concat_gdf(gdf_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preprocess #3\n",
    "- Visualise data footprints"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots(figsize=(10,10))\n",
    "temp_web = gdf_colombia.to_crs(plot_web)\n",
    "clim = np.percentile(temp_web['h_te_best_fit'].values,(2,98))\n",
    "temp_web.plot('h_te_best_fit',ax=ax,s=3,legend=True,cmap='inferno',vmin=clim[0],vmax=clim[1])\n",
    "ctx.add_basemap(ax=ax)\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### We will use the TANDEM-X Global DEM for our comparison. The resolution of the globally avaialable product is 90 m, with *horizontal* and *vertical* accuracy better than 2 to 3 m.\n",
    "- TANDEM-X DEM for the region was downloaded and preprocessed, filtered using scripts from the [tandemx](https://github.com/dshean/tandemx) repository"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dem_file = os.path.join(wd,'supporting_files/TDM1_DEM_90m_colombia_DEM_masked_aea.tif')\n",
    "hs_file = os.path.splitext(dem_file)[0]+'_hs.tif'\n",
    "dem_ds = rasterio.open(dem_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! gdaldem hillshade $dem_file $hs_file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hs_ds = rasterio.open(hs_file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gdf_on_raster(gdf,ds,ax,hs_ds=None,cmap='inferno'):\n",
    "    gdf = gdf.to_crs(ds.crs)\n",
    "    xmin,ymin,xmax,ymax = ds.bounds\n",
    "    ndv = gda_lib.get_ndv(ds)\n",
    "    img = ds.read(1)\n",
    "    img = np.ma.masked_less_equal(img,ndv)\n",
    "    clim = np.nanpercentile(img,(2,98))\n",
    "    if hs_ds:\n",
    "        hs = hs_ds.read(1)\n",
    "        ndv = gda_lib.get_ndv(hs_ds)\n",
    "        hs = np.ma.masked_less_equal(hs,ndv)\n",
    "        ax.imshow(hs,cmap='gray',extent=[xmin,xmax,ymin,ymax])\n",
    "        im = ax.imshow(img,alpha=0.6,cmap=cmap,extent=[xmin,xmax,ymin,ymax])\n",
    "        print(clim)\n",
    "    else:\n",
    "        im = ax.imshow(img,cmap=cmap,vmin=clim[0],vmax=clim[1],extent=[xmin,xmax,ymin,ymax])\n",
    "    gdf.plot('p_b',ax=ax,s=1)\n",
    "    plt.colorbar(im,ax=ax,extend='both',label='Elevation (m)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xmin,ymin,xmax,ymax = dem_ds.bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Filter points based on DEM extent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_colombia['x_atc'] = gdf_colombia['delta_time']\n",
    "gdf_colombia_dem_extent = gdf_colombia.to_crs(dem_ds.crs).cx[xmin:xmax,ymin:ymax]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots(figsize=(10,5))\n",
    "gdf_on_raster(gdf_colombia_dem_extent,dem_ds,ax,hs_ds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 1\n",
    "- This contains demonstration of elevation profile along 1 track, which has 6 beams"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Picking out 1 track\n",
    "### check with different ATLO8 inputs\n",
    "test_track = ATL08_list[3]\n",
    "print(test_track)\n",
    "test_gdf = gda_lib.ATL08_2_gdf(test_track,dataset_dict).to_crs(dem_ds.crs).cx[xmin:xmax,ymin:ymax]\n",
    "fig,ax = plt.subplots(figsize=(10,5))\n",
    "gdf_on_raster(test_gdf,dem_ds,ax,hs_ds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Working with track from 20190105 to show how we can use this to plot elevation values along the collect just by using ICESat-2\n",
    "np.unique(test_gdf['p_b'].values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Limit analysis to 1 pair beam combination\n",
    "mask = test_gdf['p_b']== '1.0_0.0'\n",
    "test_gdf_pb = test_gdf[mask]\n",
    "fig,ax = plt.subplots(figsize=(5,4))\n",
    "plot_var = test_gdf_pb['h_te_best_fit'].values\n",
    "ax.scatter(np.arange(len(plot_var)),plot_var,s=1)\n",
    "ax.set_xlabel('Along-track id')\n",
    "ax.set_ylabel('ATL08-Terrain Height')\n",
    "ax.grid('--')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#or do it for all the 6 tracks\n",
    "track_identifier = list(np.unique(test_gdf['p_b'].values))\n",
    "fig,axa = plt.subplots(3,2,figsize=(10,10))\n",
    "ax = axa.ravel()\n",
    "for idx,track in enumerate(track_identifier):\n",
    "    mask = test_gdf['p_b']== track\n",
    "    test_gdf_pb = test_gdf[mask]\n",
    "    plot_var = test_gdf_pb['h_te_best_fit'].values\n",
    "    ax[idx].scatter(np.arange(len(plot_var)),plot_var,s=1)\n",
    "    ax[idx].set_xlabel('Along-track id')\n",
    "    ax[idx].set_ylabel('ATL08-Terrain Height')\n",
    "    ax[idx].grid('--')\n",
    "    ax[idx].set_title('Track: {} Beam: {}'.format(track.split('_',15)[0],track.split('_',15)[1]))\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 2:\n",
    "- Compare ICESat-2 Elevation with that of reference DEM (in this case TANDEM-X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Sample elevations from DEM at ATLO8-locations using nearest neighbour algorithm "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del_time,elev = gda_lib.sample_near_nbor(dem_ds,gdf_colombia_dem_extent)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_colombia_dem_extent['dem_z'] = elev"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Plot elevation differences (ICESat-2 minus TANDEM-X) as a function of elevation\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_colombia_dem_extent['z_diff'] = gdf_colombia_dem_extent['h_te_best_fit'] - gdf_colombia_dem_extent['dem_z']\n",
    "fig,ax = plt.subplots(figsize=(5,4))\n",
    "# Sort elevation values\n",
    "gdf_colombia_dem_extent.sort_values(by='dem_z',inplace=True)\n",
    "gdf_colombia_dem_extent_filt = gdf_colombia_dem_extent[gdf_colombia_dem_extent['z_diff']<1000]\n",
    "ax.scatter(gdf_colombia_dem_extent_filt.dem_z.values,gdf_colombia_dem_extent_filt.z_diff.values,s=1)\n",
    "ax.set_ylim(-50,50)\n",
    "ax.set_xlabel('Elevation (TANDEM-X) (m)')\n",
    "ax.set_ylabel('Elevation difference (m)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- The difference above might be noise or real signal (\" the dates of ICESAT-2 footprints are between December to March 2018-2019, while TANDEM-X contains a mosaic of elevations between 2012-2014\")\n",
    "- It's hard to make out anything from the above plot, let's try a box plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dem_bins = list(np.arange(0,5500,500))\n",
    "# mask out differences larger than 100 m ?\n",
    "filt_lim = (-100,100)\n",
    "mask = (gdf_colombia_dem_extent['z_diff']<=100) & (gdf_colombia_dem_extent['z_diff']>=-100)\n",
    "gdf_colombia_dem_extent_filt_box = gdf_colombia_dem_extent[mask]\n",
    "gdf_colombia_dem_extent_filt_box['bin'] = pd.cut(gdf_colombia_dem_extent_filt_box['dem_z'],bins=dem_bins)\n",
    "fig,ax = plt.subplots(figsize=(5,4))\n",
    "gdf_colombia_dem_extent_filt_box.boxplot(column='z_diff',by='bin',ax=ax)\n",
    "ax.set_xlabel('Elevation (TANDEM-X) (m)')\n",
    "ax.set_xticklabels(dem_bins)\n",
    "#ax.set_ylabel('Elevation difference (m)')\n",
    "ax.set_title('')\n",
    "ax.set_ylabel('ICESat-2 minus TANDEM-X (m)')\n",
    "#plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- The x labels in the plot are lower intervals of boxplots, we see that the median differences are close to zero for most elevation ranges with a maximum difference of -10 m. Also, we see a constant negative bias in all the elevation difference. This might be due to a bias present between the 2 sources. This bias maybe due to offset between the 2 datasets which might come down after coregistration."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 3\n",
    "- Application of ICESat-2 as control surface for DEMs coregistration\n",
    "- Or, to find offsets and align ICESat-2 tracks to a control surface"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Going fancy, include only if you want to :)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Application of ICESat-2 as control for DEM co-registration ?\n",
    "- Can use point cloud alignment techniques to align DEMs to points, for now as a starting point we can use the transformation matrix to inform on the horizontal and vertical offset between ICESat-2 tracks and DEMs\n",
    "- We will be using a flavor of Iterative Closest Point alignment algorithm, implemented in [Ames Stereo Pipeline](https://github.com/NeoGeographyToolkit/StereoPipeline)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_colombia_dem_extent.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Save the geodataframe in the specified way as expected by Ames Stereo Pipline\n",
    "icesat2_pc = '/home/jovyan/icesat2/icesat2_colombia_pc.csv' \n",
    "gdf_colombia_dem_extent[['latitude','longitude','h_te_best_fit']].to_csv(icesat2_pc,header=False,index=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_colombia_dem_extent.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Save the geodataframe in the specified way as expected by Ames Stereo Pipline\n",
    "icesat2_pc = '/home/jovyan/icesat2/icesat2_colombia_pc.csv'\n",
    "pc_rename_dict = {'latitude':'lat','longitude':'lon','h_te_best_fit':'height_above_datum'}\n",
    "gdf_colombia_dem_extent = gdf_colombia_dem_extent.rename(columns=pc_rename_dict)\n",
    "#gdf_colombia_dem_extent['height_above_datum'] = gdf_colombia_dem_extent['h_te_best_fit']\n",
    "\n",
    "gdf_colombia_dem_extent[['lon','lat','height_above_datum']].to_csv(icesat2_pc,header=True,index=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!exportPATH=\"/home/jovyan/icesat2/StereoPipeline/bin:$PATH\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "! ls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "align_fol = '/home/jovyan/icesat2/align/run'\n",
    "#max-displacement is set to 10, given ICESat-2 reported operational accuracy\n",
    "pc_align_opts=\"--csv-format '1:lon 2:lat 3:height_above_datum' --max-displacement 10 --save-transformed-source-points --alignment-method point-to-point --datum WGS84\"\n",
    "!/home/jovyan/icesat2/StereoPipeline/bin/pc_align $pc_align_opts $icesat2_pc $dem_file -o $align_fol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Alignment results suggest that there is an offset of ~5.4 m between the ICESat-2 points and TANDEM-X DEM, so that could have contributed to the offsets which we see above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##Lets rerun the analysis with the new DEM to see if the alignment improved anything or not\n",
    "## Regrid the transformed pointcloud into DEM at 90 m posting\n",
    "!/home/jovyan/icesat2/StereoPipeline/bin/point2dem --tr 90 --t_srs EPSG:32618 $align_fol-trans_source.tif"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_colombia_dem_extent = gdf_colombia_dem_extent.loc[:,~gdf_colombia_dem_extent.columns.duplicated()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf_colombia_dem_extent['height_above_datum'].values[5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trans_dem_file = '/home/jovyan/icesat2/align/run-trans_source-DEM.tif'\n",
    "trans_dem_ds = rasterio.open(trans_dem_file)\n",
    "del_time,elev = gda_lib.sample_near_nbor(trans_dem_ds,gdf_colombia_dem_extent)\n",
    "gdf_colombia_dem_extent['trans_dem_z'] = elev\n",
    "dem_bins = list(np.arange(0,5500,500))\n",
    "# mask out differences larger than 100 m ?\n",
    "filt_lim = (-100,100)\n",
    "gdf_colombia_dem_extent['trans_z_diff'] = gdf_colombia_dem_extent.height_above_datum - gdf_colombia_dem_extent.trans_dem_z\n",
    "\n",
    "mask = (gdf_colombia_dem_extent['trans_z_diff']<=100) & (gdf_colombia_dem_extent['trans_z_diff']>=-100)\n",
    "gdf_colombia_dem_extent_filt_box = gdf_colombia_dem_extent[mask]\n",
    "gdf_colombia_dem_extent_filt_box['bin'] = pd.cut(gdf_colombia_dem_extent_filt_box['dem_z'],bins=dem_bins)\n",
    "fig,ax = plt.subplots(figsize=(5,4))\n",
    "gdf_colombia_dem_extent_filt_box.boxplot(column='trans_z_diff',by='bin',ax=ax)\n",
    "ax.set_xlabel('Elevation (TANDEM-X) (m)')\n",
    "ax.set_xticklabels(dem_bins)\n",
    "ax.set_title('')\n",
    "ax.set_ylabel('ICESat-2 minus TANDEM-X DEM after coregistration (m)')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- We see that after coregistration, the bias reduces to an extent. Note that this is a very preliminary analysis, results will be better after filtering the ATL08 points based on quality metrics and finding truly static surfaces (snow free during acquisiton time of ICESat-2 points)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}