access_atl08.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cd7a46d3",
   "metadata": {},
   "source": [
    "# S3 Access ICESat-2 ATL08 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4ed75e2b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# hide warning messages \n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "# import python modules\n",
    "import geopandas as gpd\n",
    "import requests\n",
    "import pandas as pd\n",
    "import datetime as dt \n",
    "import h5py\n",
    "import numpy as np\n",
    "import s3fs\n",
    "from os import path, remove\n",
    "from posixpath import splitext\n",
    "from requests.adapters import HTTPAdapter, Retry\n",
    "import concurrent.futures"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f66c9eb",
   "metadata": {},
   "source": [
    "The following input parameters are needed, passed through the front-end."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "2b6ceb2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# start time\n",
    "start_date = dt.datetime(2017, 1, 1) \n",
    "# end time\n",
    "end_date = dt.datetime(2020, 1, 31)\n",
    "# atl08 variables of interests as a list\n",
    "variables = ['land_segments/canopy/h_canopy', 'land_segments/canopy/canopy_h_metrics',\n",
    "             'signal_photons/ph_h']\n",
    "# path to polygon geojson file define subset area\n",
    "poly_f = 'poly.json'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1532fce",
   "metadata": {},
   "source": [
    "Let's define some additional variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "52fa7e44",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ICESat-2 reference orbits\n",
    "orbits = gpd.read_file('orbits.json')\n",
    "# ATL08 regions/granule boundaries\n",
    "region = gpd.read_file('region.json')\n",
    "# nsidc s3 credentials\n",
    "nsidc_s3 = f\"https://data.nsidc.earthdatacloud.nasa.gov/s3credentials\"\n",
    "# cmr api url\n",
    "cmrurl='https://cmr.earthdata.nasa.gov/search/' \n",
    "# ATL08 (Earthdata Cloud) concept_id\n",
    "concept_id = 'C2153574670-NSIDC_CPRD' # 'C2144424132-NSIDC_ECS' is the public concept id\n",
    "# CMR datetime format\n",
    "dt_format = '%Y-%m-%dT%H:%M:%SZ'\n",
    "# six beams of ICESat-2\n",
    "beams = ['gt1l', 'gt1r', 'gt2l', 'gt2r', 'gt3l', 'gt3r']\n",
    "# core variables to include in every subset request\n",
    "headers = ['land_segments/latitude', 'land_segments/longitude', \n",
    "           'land_segments/segment_id_beg', 'land_segments/segment_id_end',\n",
    "           'beam', 'granule']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8807a1e",
   "metadata": {},
   "source": [
    "Let's create a session to avoid 500 request errors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "f08ad763",
   "metadata": {},
   "outputs": [],
   "source": [
    "# creating a backoff\n",
    "s = requests.Session()\n",
    "retries = Retry(total=3, backoff_factor=0.1, status_forcelist=[ 500, 502, 503, 504 ])\n",
    "s.mount('https://', HTTPAdapter(max_retries=retries))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a47f8680",
   "metadata": {},
   "source": [
    "Let's find the intersecting ICESat-2 reference orbits intersecting with the polygon (`poly`). For intersection to work correctly, geopandas requires [pygeos](https://pygeos.readthedocs.io/) module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "d0ed4a40",
   "metadata": {},
   "outputs": [],
   "source": [
    "poly = gpd.read_file(poly_f)\n",
    "poly['geometry'] = poly.geometry.buffer(0.3, cap_style=3)\n",
    "inter = gpd.sjoin(orbits, poly, how='inner', op='intersects')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "129805d4",
   "metadata": {},
   "source": [
    "The orbit column in the `inter` dataframe will have a four-digit reference ground track number of the intersecting orbits. We will now search CMR for ATL08 granules by passing this number to a CMR API request. This is necessary since the CMR bounding geometries of ATL08 are not representative of actual data bounds.\n",
    "\n",
    "Accessing s3 links of atl08 via CMR requires authentication as it is not yet publicly available. Create a NASA EDL token first by going to `https://urs.earthdata.nasa.gov/users/<username>/user_tokens` and generate a token if you don't have one already. Save the token to the `~/.cmr_token` file in your home directory.\n",
    "\n",
    "Now, let's define functions to search for ATL08 granules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b45b87d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# function to search CMR \n",
    "def granule_search(concept_id: str, temporal: str, granule_id:str):\n",
    "    page_num = 1\n",
    "    granule_arr = []\n",
    "    while True:\n",
    "         # defining parameters\n",
    "        cmr_param = {\n",
    "            \"collection_concept_id\": concept_id, \n",
    "            \"temporal\": temporal,\n",
    "            \"page_size\": 2000,\n",
    "            \"page_num\": page_num,\n",
    "            \"options[producerGranuleId][pattern]\": \"true\",\n",
    "            \"producerGranuleId\": granule_id        }\n",
    "\n",
    "        granulesearch = cmrurl + 'granules.json'\n",
    "        token_path = path.join(path.expanduser(\"~\"), '.cmr_token')\n",
    "        response = requests.get(granulesearch, params=cmr_param)\n",
    "        granules = response.json()['feed']['entry']\n",
    "        if granules:\n",
    "            for g in granules:\n",
    "                # Get URL of HDF5 files\n",
    "                for links in g['links']:\n",
    "                    if links['href'].startswith('s3://') and links['href'].endswith('.h5'):\n",
    "                        granule_arr.append(links['href'])\n",
    "\n",
    "            page_num += 1\n",
    "        else: \n",
    "            break\n",
    "            \n",
    "    return granule_arr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50290061",
   "metadata": {},
   "source": [
    "We now pass time bounds and four-digit ICESat-2 reference orbit numbers to the CMR search function defined above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4a436a8c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# CMR time bounds\n",
    "temporal = start_date.strftime(dt_format) + ',' + end_date.strftime(dt_format)\n",
    "alt08_f = []\n",
    "for o_n in inter.orbit:\n",
    "    alt08_f.extend(granule_search(concept_id, temporal, f\"*_{str(o_n).zfill(4)}*\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e43caf85",
   "metadata": {},
   "source": [
    "The list `alt08_f` will now have S3 links to ATL08 granules.\n",
    "\n",
    "Each ATL08 granules covers about 1/14th of an orbit, called \"region\", identified by a two-digit number. Each region is defined by latitude boundaries and whether the orbit is ascending or descending. Let's find all the regions that intersects the polygon `poly`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c9d2b268",
   "metadata": {},
   "outputs": [],
   "source": [
    "regions = []\n",
    "for r_inter in gpd.sjoin(region, poly, how='inner', op='intersects').region:\n",
    "    for r in r_inter.split(','):\n",
    "        regions.append(r)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ebec4a3",
   "metadata": {},
   "source": [
    "Now, we will only keep the granules that belongs to the regions identified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "52d71f4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "alt08_files = []\n",
    "for s3_url in alt08_f:\n",
    "    granule_region = path.basename(s3_url).split(\"_\")[2][-2:]\n",
    "    if granule_region in regions:\n",
    "        alt08_files.append(s3_url)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3cf896ae",
   "metadata": {},
   "source": [
    "We will now get S3 temporary credentials and define the S3FS S3FileSystem object. Make sure that `.netrc` file is defined as described here: https://wiki.earthdata.nasa.gov/display/EL/How+To+Access+Data+With+cURL+And+Wget."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "027e2a39",
   "metadata": {},
   "outputs": [],
   "source": [
    "# get s3 credentials\n",
    "s3credentials = s.get(nsidc_s3).json()\n",
    "# defining S3FS object\n",
    "fs_s3 = s3fs.S3FileSystem(anon=False, \n",
    "                          key=s3credentials['accessKeyId'], \n",
    "                          secret=s3credentials['secretAccessKey'], \n",
    "                          token=s3credentials['sessionToken'])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4dcd9088",
   "metadata": {},
   "source": [
    "We will now define a function to spatially subset the ATL08 granules. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "2311748d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_ph_i(beg, end, ph_segment_id, hf_var):\n",
    "    z1 = np.where(ph_segment_id==beg)[0][0]\n",
    "    z2 = np.where(ph_segment_id==end)[0][-1]\n",
    "    return (hf_var[z1:z2+1].tolist())\n",
    "\n",
    "def subset_atl08(t, var_l, poly, s3_url):\n",
    "    photon = False\n",
    "    if any(v.startswith('signal_photons') for v in var_l):\n",
    "        photon = True\n",
    "        \n",
    "    with fs_s3.open(s3_url, mode='rb') as fh:\n",
    "        with h5py.File(fh) as hf:\n",
    "            granule_f = path.basename(s3_url)\n",
    "            out_csv = f'{splitext(granule_f)[0]}_{int(dt.datetime.utcnow().timestamp())}.csv'\n",
    "            for var in list(hf.keys()):\n",
    "                if var.startswith('gt'):\n",
    "                    if 'land_segments' in list(hf[var].keys()):\n",
    "                        lat = hf[f'{var}/{var_l[0]}'][:]\n",
    "                        lon = hf[f'{var}/{var_l[1]}'][:]\n",
    "                        df = pd.DataFrame({var_l[0]: lat, var_l[1]: lon})\n",
    "                        if photon:\n",
    "                            df[var_l[2]] = hf[f'{var}/{var_l[2]}'][:]\n",
    "                            df[var_l[3]] = hf[f'{var}/{var_l[3]}'][:]\n",
    "                        gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df[var_l[1]], df[var_l[0]])) \n",
    "                        gdf_poly = gdf[gdf['geometry'].within(poly.geometry[0])]   \n",
    "                        if not gdf_poly.empty:\n",
    "                            if not path.exists(out_csv):\n",
    "                                with open(out_csv, \"w\") as f:\n",
    "                                    f.write(','.join(var_l)+'\\n')\n",
    "                        \n",
    "                            gdf_poly['beam'] = np.repeat(str(var), len(gdf_poly.index)).tolist()\n",
    "                            gdf_poly['granule'] = np.repeat(str(granule_f), len(gdf_poly.index)).tolist()\n",
    "                            \n",
    "                            if photon:\n",
    "                                ph_segment_id = hf[f'{var}/signal_photons/ph_segment_id'][:]\n",
    "\n",
    "                            for v in var_l[t:]:\n",
    "                                if v not in var_l[:t]:\n",
    "                                    gdf_poly[v] = None\n",
    "\n",
    "                            # retrieving variables of interests\n",
    "                            for _, df_gr in gdf_poly.groupby((gdf_poly.index.to_series().diff() > 1).cumsum()):\n",
    "                                i = df_gr.index.min()\n",
    "                                j = df_gr.index.max()\n",
    "                                for v in var_l[t:]:\n",
    "                                    if v not in var_l[:t]:\n",
    "                                        if v.startswith('signal_photons'):\n",
    "                                            gdf_poly.loc[i:j, (v)][:] = [find_ph_i(x1, x2, \n",
    "                                                                                   ph_segment_id, \n",
    "                                                                                   hf[f'{var}/{v}']) \n",
    "                                                                         for x1, x2 in zip(\n",
    "                                                    gdf_poly[var_l[2]], gdf_poly[var_l[3])]\n",
    "                                        else:\n",
    "                                            gdf_poly.loc[i:j, (v)][:] = hf[f'{var}/{v}'][i:j+1].tolist()\n",
    "                                    \n",
    "\n",
    "                            # saving the output file\n",
    "                            gdf_poly.to_csv(out_csv, mode='a', index=False, header=False, columns=var_l)\n",
    "    if path.exists(out_csv):                         \n",
    "        return out_csv"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6fa4d3e",
   "metadata": {},
   "source": [
    "To the above function `subset_atl08`, we will pass a list `headers` and `variables` containing the variables of interest, a geopandas dataframe `poly_gpd` containing subset area polygon, and a list `alt08_files` containing the list of S3 links to ATL08 files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3ebbaafc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read the polygon json defining the subset area\n",
    "poly_gpd = gpd.read_file(poly_f)\n",
    "# subset using asynchronous calls\n",
    "futures = []\n",
    "subset_f = []\n",
    "subset_pool = concurrent.futures.ThreadPoolExecutor()\n",
    "for s3_url in sorted(alt08_files):\n",
    "    future = subset_pool.submit(subset_atl08, len(headers), headers+variables, poly_gpd, s3_url)\n",
    "    futures.append(future)\n",
    "futures, _ = concurrent.futures.wait(futures)\n",
    "for future in futures:\n",
    "    result = future.result()\n",
    "    subset_f.append(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7f13c51",
   "metadata": {},
   "source": [
    "The variable `subset_f` will now contain a list of CSV files that were created in the above step.  These CSV files contains the subset data. Now, we will read these subset CSV files and print out the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b89cc93e",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# reading the CSV\n",
    "subset_f = [x for x in subset_f if x is not None]\n",
    "subset_df = pd.concat(pd.read_csv(f, header=0) for f in subset_f)\n",
    "subset_df.reset_index(inplace=True)\n",
    "# deleting the CSV\n",
    "for f in subset_f:\n",
    "    if path.isfile(f):\n",
    "        remove(f)\n",
    "# printing the pandas dataframe as json\n",
    "subset_df.to_json()"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}