geemap_ICESat2_OA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "deaa33ed-ff93-4326-ae72-fd69fcf374b1",
   "metadata": {},
   "source": [
    "# ICESat-2 Meets Google Earth Engine Map\n",
    "\n",
    "If you run this in a new computational environment for the first time \n",
    "(including on Binder, Colab, etc.) you need to authenticate google earth engine. \n",
    "Open a terminal and run\n",
    "~~~\n",
    "$earthengine authenticate\n",
    "~~~\n",
    "Then follow the instructions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4855f3e7-3283-43ff-9339-2a9f9fe09f06",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pylab as plt\n",
    "import numpy as np\n",
    "import ee\n",
    "import geemap\n",
    "import datetime\n",
    "import holoviews as hv\n",
    "from holoviews import opts\n",
    "hv.extension('bokeh', 'matplotlib')\n",
    "\n",
    "from oautils import dataCollector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26edc99a-9e85-49cd-a691-12b6f4fd26b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # uncomment if you need to authenticate Google Earth Engine\n",
    "ee.Authenticate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57e55864-e9fc-45a0-9f87-a5d8b8e447ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b65219e6-faeb-41e7-8897-f2d29616eb73",
   "metadata": {},
   "source": [
    "### Specify the parameters for the ICESat-2 data\n",
    "Browse OpenAltimetry to find a track / region of interest! Here's a few examples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b211816-3e8a-418a-a7ba-1ee19fbca1b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "##################################################################\n",
    "# amery icesat-2 doline by grounding line\n",
    "# kwargs = {'latlims': [-71.3364, -71.280],\n",
    "#           'lonlims': [71.450, 71.536],\n",
    "#           'track': 721,\n",
    "#           'beam': 'gt3r'}\n",
    "# date = '2019-11-13'\n",
    "# date = '2019-05-15'\n",
    "\n",
    "##################################################################\n",
    "# # weird amery 'lake' / blue ice? \n",
    "# kwargs = {'latlims': [-72.44, -72.418],\n",
    "#           'lonlims': [67.40, 67.44],\n",
    "#           'track': 81,\n",
    "#           'beam': 'gt3r'} # also in gt3l (strong beam)\n",
    "# date = '2019-01-02'\n",
    "\n",
    "##################################################################\n",
    "# lake on Nivlisen ice shelf\n",
    "# http://openaltimetry.org/data/api/icesat2/atl03?date=2020-01-16&minx=12.107692195781404&miny=-70.34956862465471&maxx=12.426364789894341&maxy=-70.2449105354736&trackId=312&beamName=gt3r&beamName=gt3l&beamName=gt2r&beamName=gt2l&beamName=gt1r&beamName=gt1l&outputFormat=json\n",
    "kwargs = {'latlims': [-70.34957, -70.24491],\n",
    "          'lonlims': [12.10769, 12.42636],\n",
    "          'track': 312,\n",
    "          'beam': 'gt2r'} # also in gt2l (weak beam)\n",
    "date = '2020-01-16'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c93760bf-d496-40e3-9d84-b68732f80899",
   "metadata": {},
   "source": [
    "### Now download the data\n",
    "This uses the dataCollector class defined in oautils.py. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b57e6e0-63bb-46bb-afc7-d5e9824c2560",
   "metadata": {},
   "outputs": [],
   "source": [
    "is2data = dataCollector(date=date,**kwargs)\n",
    "is2data.requestData()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e57df456-bf5a-4fc9-add8-ade7c2d7c774",
   "metadata": {},
   "source": [
    "### The ATL03 / ATL06 / ATL08 data are now available as pandas dataframes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e97df09-1925-4be1-bd74-c5254dcda6b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "is2data.atl06"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9cffb997-9563-40ed-adb1-b8f2a024333a",
   "metadata": {},
   "source": [
    "### Plot the ATL03 and ATL06/ATL08 data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf67e481-eb7c-4629-b110-353dbcd7d1cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "confdict = {'Noise': -1.0, 'Buffer': 0.0, 'Low': 1.0, 'Medium': 2.0, 'High': 3.0}\n",
    "is2data.atl03['conf_num'] = [confdict[x] for x in is2data.atl03.conf]\n",
    "is2data.atl08['canopy_h'] = is2data.atl08.h + is2data.atl08.canopy\n",
    "atl03scat = hv.Scatter(is2data.atl03, 'lat', vdims=['h', 'conf_num'], label='ATL03')\\\n",
    "            .opts(color='conf_num', alpha=1, cmap='dimgray_r')\n",
    "atl06line = hv.Curve(is2data.atl06, 'lat', 'h', label='ATL06')\\\n",
    "            .opts(color='r', alpha=0.5, line_width=3)\n",
    "atl08line = hv.Curve(is2data.atl08, 'lat', 'h', label='ATL08')\\\n",
    "            .opts(color='b', alpha=1, line_width=1)\n",
    "# atl08scat = hv.Scatter(is2data.atl08, 'lat', 'canopy_h', label='ATL08 Canopy')\n",
    "# atl08scat = atl08scat.opts(alpha=1, color='b', size=4)\n",
    "hrange = is2data.atl06.h.max() - is2data.atl06.h.min()\n",
    "overlay = (atl03scat * atl06line * atl08line).opts(\n",
    "    height=500, \n",
    "    width=900,\n",
    "    xlabel='latitude', \n",
    "    ylabel='elevation', \n",
    "    title='ICESat-2 track %d %s on %s' % (kwargs['track'],kwargs['beam'].upper(),date),\n",
    "    legend_position='bottom_right',\n",
    "    ylim=(is2data.atl06.h.min()-hrange ,is2data.atl06.h.max()+hrange),\n",
    "    xlim=(is2data.atl06.lat.min(), is2data.atl06.lat.max())\n",
    ")\n",
    "overlay"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cfa64464-0e6c-4769-9ed3-679a25338c6a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# fig = plt.figure(figsize=[9,4.5],dpi=100)\n",
    "# ax = fig.add_subplot(1,1,1)\n",
    "# ax.scatter(is2data.atl03.lat,is2data.atl03.h,s=10,c='k',alpha=0.1,edgecolors='none')\n",
    "# h2, = ax.plot(is2data.atl06.lat,is2data.atl06.h,c='r',lw=3, alpha=0.5, label='ATL06 land ice elevation')\n",
    "# h3, = ax.plot(is2data.atl08.lat,is2data.atl08.h,c='c',lw=0.8, ls='--', alpha=1, label='ATL08 land height')\n",
    "# h4 = ax.scatter(is2data.atl08.lat,is2data.atl08.canopy,s=15,c='c', alpha=1, label='ATL08 canopy')\n",
    "# h1 = ax.scatter([],[],s=20,c='k',alpha=0.5,edgecolors='none',label='ATL03 photon cloud')\n",
    "# # h2, = ax.plot([-9999, -9998],[-9999, -9998],c='r',lw=1.3, label='ATL06 land ice elevation')\n",
    "# ax.set_xlim(kwargs['latlims'])\n",
    "# hrange = is2data.atl06.h.max() - is2data.atl06.h.min()\n",
    "# ax.set_ylim(is2data.atl06.h.min() - 0.5*hrange ,is2data.atl06.h.max() + 0.5*hrange)\n",
    "# ax.set_xlabel('latitude',fontsize=8)\n",
    "# ax.set_ylabel('elevation [m]',fontsize=8)\n",
    "# ax.tick_params(axis='both', labelsize=7)\n",
    "# ax.legend(handles=[h1,h2,h3,h4], loc='lower right', fontsize=7)\n",
    "# ax.set_title('ICESat-2 track %d %s on %s' % (kwargs['track'],kwargs['beam'].upper(),date),fontsize=9)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bcb412d-bd86-4da7-8c70-ab81f463a2c3",
   "metadata": {},
   "source": [
    "### How long is this ground track?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bcf12521-228f-450c-819d-6b6dc1a93bfa",
   "metadata": {},
   "outputs": [],
   "source": [
    "lat1, lat2 = is2data.atl06.lat[0], is2data.atl06.lat.iloc[-1]\n",
    "lon1, lon2 = is2data.atl06.lon[0], is2data.atl06.lon.iloc[-1]\n",
    "\n",
    "def dist_latlon2meters(lat1, lon1, lat2, lon2):\n",
    "    # returns the distance between two lat/lon coordinate points along the earth's surface in meters\n",
    "    R = 6371000\n",
    "    def deg2rad(deg):\n",
    "        return deg * (np.pi/180)\n",
    "    dlat = deg2rad(lat2-lat1)\n",
    "    dlon = deg2rad(lon2-lon1)\n",
    "    a = np.sin(dlat/2) * np.sin(dlat/2) + np.cos(deg2rad(lat1)) * np.cos(deg2rad(lat2)) * np.sin(dlon/2) * np.sin(dlon/2)\n",
    "    c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1-a))\n",
    "    return R * c\n",
    "\n",
    "dist_xatc = dist_latlon2meters(lat1, lon1, lat2, lon2)\n",
    "print('The ground track is %d meters long.' % np.round(dist_xatc))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00b87af8-f609-4180-87e4-a84d2c102ffd",
   "metadata": {},
   "source": [
    "### Add the ICESat-2 ground track to the map\n",
    "We can use google earth satellite imagery as the basemap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa33c611-e2d5-47cf-ba9b-f9744662c3f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "gtx_geom = ee.Geometry.LineString(coords=list(zip(is2data.atl06.lon,is2data.atl06.lat)), proj='EPSG:4326', geodesic=True)\n",
    "Map = geemap.Map(center=(40, -100), zoom=4)\n",
    "Map.add_basemap('HYBRID')\n",
    "Map.addLayer(ee.FeatureCollection(gtx_geom),{'color': 'red'},'ground track')\n",
    "center_lat = (lat1 + lat2) / 2\n",
    "center_lon = (lon1 + lon2) / 2\n",
    "Map.setCenter(center_lon, center_lat, 7);\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f515d4e-9111-4fb0-8aab-0fb364029270",
   "metadata": {},
   "source": [
    "### Let's look for Sentinel-2 and Landsat 8 images!\n",
    "These live in image collections on Google Earth Engine. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad33783d-28c9-4c52-ad3f-971d1ff1f526",
   "metadata": {},
   "outputs": [],
   "source": [
    "collection_name1 = 'COPERNICUS/S2_SR'  # Sentinel-2 earth engine collection \n",
    "# https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2_SR\n",
    "\n",
    "collection_name2 = 'LANDSAT/LC08/C01/T2'  # Landsat 8 earth engine collection \n",
    "# https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C01_T2\n",
    "# Note: Landsat 8 ingestion into Earth Engine seems to not have reached Antarctica yet, so using raw scenes...\n",
    "\n",
    "# Specify for how many days around the ICESat-2 overpass we want to look \n",
    "days_buffer_imagery = 15"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4c2fa38-460d-4810-9a8f-d5cd546924dc",
   "metadata": {},
   "source": [
    "### Query for images that overlap spatially with - and were acquired close in time to - the ICESat-2 overpass."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6dc541ec-bd2a-4b09-bd31-b25155ff72f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "date_requested = date\n",
    "dateformat = '%Y-%m-%d'\n",
    "datetime_requested = datetime.datetime.strptime(date_requested, dateformat)\n",
    "search_start = (datetime_requested - datetime.timedelta(days=days_buffer_imagery)).strftime(dateformat)\n",
    "search_end = (datetime_requested + datetime.timedelta(days=days_buffer_imagery)).strftime(dateformat)\n",
    "\n",
    "# the point of interest (center of the track) as an Earth Engine Geometry\n",
    "poi = ee.Geometry.Point(center_lon, center_lat)\n",
    "\n",
    "# the region of interest around it\n",
    "buffer_around_center_meters = dist_xatc*0.8\n",
    "roi = poi.buffer(buffer_around_center_meters)\n",
    "\n",
    "# the collection to query: \n",
    "# 1) merge Landsat 8 and Sentinel-2 collections\n",
    "# 2) filter by acquisition date\n",
    "# 3) filter by the point of interest\n",
    "# 4) sort by acquisition date\n",
    "collection = ee.ImageCollection(collection_name1) \\\n",
    "    .merge(ee.ImageCollection(collection_name2)) \\\n",
    "    .filterDate(search_start, search_end) \\\n",
    "    .filterBounds(poi) \\\n",
    "    .sort('system:time_start') \n",
    "\n",
    "# clip to the region of interest\n",
    "def clip(img): return img.clip(roi)\n",
    "collection = collection.map(clip);\n",
    "\n",
    "# get info \n",
    "n_imgs = collection.size().getInfo()\n",
    "info = collection.getInfo()\n",
    "print('number of images +/- %d days of date requested: %d' % (days_buffer_imagery, n_imgs))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "981d8c82-7d6f-4e40-815f-139e7ea7e8ff",
   "metadata": {},
   "source": [
    "### Add all the images that were found to the Map as individual layers\n",
    "If the image is Landsat 8, this pan-sharpens the RGB image to give better spatial resolution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aeb0467a-4d46-4c0b-a239-1df9ed6b84a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "Map = geemap.Map(center=(40, -100), zoom=4)\n",
    "Map.add_basemap('HYBRID')\n",
    "\n",
    "# get a list of images to access them in sequence (there may be a better way to do this...)\n",
    "list_of_images = collection.toList(n_imgs)\n",
    "for i in range(n_imgs):\n",
    "    \n",
    "    # get the relevant info\n",
    "    thisDate = datetime.datetime.fromtimestamp(info['features'][i]['properties']['system:time_start']/1e3)\n",
    "    dtstr = thisDate.strftime(dateformat)\n",
    "    dt = (thisDate - datetime_requested).days\n",
    "    ID = info['features'][i]['id']\n",
    "    rel = 'before' if dt<0 else 'after'\n",
    "    print('%02d: %s (%3d days %s ICESat-2): %s' % (i, dtstr, np.abs(dt), rel, ID))\n",
    "    \n",
    "    # normalize the image to [0.0, 1.0] range, based on max and min in the rgb bands\n",
    "    img = ee.Image(list_of_images.get(i))\n",
    "    rgb = img.select('B4', 'B3', 'B2')\n",
    "    rgbmax = rgb.reduce(ee.Reducer.max()).reduceRegion(reducer=ee.Reducer.max(), geometry=rgb.geometry(), bestEffort=True, maxPixels=1e6)\n",
    "    rgbmin = rgb.reduce(ee.Reducer.min()).reduceRegion(reducer=ee.Reducer.min(), geometry=rgb.geometry(), bestEffort=True, maxPixels=1e6)\n",
    "    rgb = rgb.unitScale(ee.Number(rgbmin.get('min')), ee.Number(rgbmax.get('max'))).clamp(0.0, 1.0)\n",
    "    \n",
    "    # only show the first layer, then can toggle others on in map\n",
    "    show_layer = True if i==0 else False\n",
    "    \n",
    "    # if the image is Landsat 8, then pan-sharpen the image\n",
    "    if 'LANDSAT' in ID: \n",
    "        pan = img.select('B8').unitScale(ee.Number(rgbmin.get('min')), ee.Number(rgbmax.get('max'))).clamp(0.0, 1.0)\n",
    "        huesat = rgb.rgbToHsv().select('hue', 'saturation')\n",
    "        upres = ee.Image.cat(huesat, pan).hsvToRgb()\n",
    "        Map.addLayer(upres, {'min': 0, 'max': 1}, name='%02d: %s (L8)'%(i,dtstr), shown=show_layer)\n",
    "    \n",
    "    # if the image is Sentinel-2, then just plot the\n",
    "    else:\n",
    "        Map.addLayer(rgb, {'min': 0, 'max': 1}, name='%02d: %s (S2)'%(i,dtstr), shown=show_layer)\n",
    "\n",
    "# show the ground track and center the map on it\n",
    "Map.addLayer(ee.FeatureCollection(gtx_geom),{'color': 'red'},'ground track')\n",
    "Map.centerObject(gtx_geom)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2cef8c6-3bb8-49a0-af61-20f651023b25",
   "metadata": {},
   "source": [
    "## Finally, display the imagery and the ground track on the map. \n",
    "Click on the Layer button on the top right to toggle the visibility of the images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c442f167-8ab3-4d53-94b8-1968daaa3f7f",
   "metadata": {},
   "outputs": [],
   "source": [
    "Map"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ee-map-hackweek",
   "language": "python",
   "name": "ee-map-hackweek"
  },
  "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.9.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}