python/how-tos/how-to-find-and-access-GEDI-data_earthaccess.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f5c09f31-9ee1-43f5-9921-bdae6d5172b8",
   "metadata": {},
   "source": [
    "# How to: Find and Access GEDI  Data\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4205b684-b1f3-4674-a27c-4eb191c6154c",
   "metadata": {},
   "source": [
    "This notebook will provide a guide on how to find and access L1 and L2 Global Ecosystem Dynamics Investigation (GEDI) V2 data providing high-resolution laser ranging of Earth’s forests and topography from the International Space Station (ISS). Currently, there are two methods of finding and accessing L1 and L2  GEDI Version 2 data products:\n",
    "\n",
    "1. [NASA's Earthdata Search](https://search.earthdata.nasa.gov/search)\n",
    "2. [NASA's CMR API](https://cmr.earthdata.nasa.gov/search/site/docs/search/api.html)\n",
    "\n",
    "This notebook will explain how to find and access L1 and L2 GEDI V2 data stored in Earthdata Cloud programmatically. While you can directly work with NASA's CMR API, [`earthaccess`](https://github.com/nsidc/earthaccess ) python library, mentioned in this notebook, provides a user-friendly way to login, search, and download or stream NASA Earth science data available in [Common Metadata Repository (CMR)](https://www.earthdata.nasa.gov/eosdis/science-system-description/eosdis-components/cmr). "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ecf2e23d-d60e-4bd4-9c32-bd9ee53979e1",
   "metadata": {},
   "source": [
    "Let's start with importing the required packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54a3417a-4017-42f9-ae77-5eb328d95bb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, json\n",
    "import earthaccess\n",
    "import pandas as pd\n",
    "import geopandas as gp\n",
    "import h5py\n",
    "from pprint import pprint\n",
    "from shapely import Polygon\n",
    "from shapely.geometry import Point\n",
    "from shapely.geometry.polygon import orient\n",
    "from datetime import datetime\n",
    "\n",
    "os.chdir('../..')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ef20e1d-07d2-4a78-88df-530259af28b9",
   "metadata": {},
   "source": [
    "## Authentication"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2abb3f8-b34d-4e1f-af5a-589d819e138e",
   "metadata": {},
   "source": [
    "To access or download NASA Earth data, you need a .netrc file containing your NASA Earthdata Login information is needed. You can manually create a .netrc file but you can use earthaccess package for easier authentication.  `earthaccess.login()` function is used to authenticate with NASA Earthdata Login credentials stored in a .netrc file. This function will prompt you to enter your NASA Earthdata username and password to create the .netrc file if it doesn't already exist, and then uses your account information for authentication. \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "292c484b-3682-48f3-a145-5b5a57065f19",
   "metadata": {},
   "outputs": [],
   "source": [
    "earthaccess.login(persist=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd4805c8-84ac-4675-b7b4-75e77b01a55c",
   "metadata": {},
   "source": [
    "## Search for GEDI Collections\n",
    "\n",
    "GEDI level 1 & level 2 data products are hosted by the Land Processes Distributed Active Archive Center (LP DAAC), while GEDI L3 & L4 are distributed by Oak Ridge National Laboratory (ORNL DAAC). In this example, we will use the cloud-hosted [GEDI L2B Canopy Cover and Vertical Profile Metrics Data Global Footprint Level (GEDI02_B)](https://lpdaac.usgs.gov/products/gedi02_bv002/) to find data, but the same routine can be used to access other products. To find the data we will use the `earthaccess` Python library. `earthaccess` searches [NASA's Common Metadata Repository (CMR)](https://cmr.earthdata.nasa.gov/search), a metadata system that catalogs Earth Science data and associated metadata records. `collection_query` from `earthaccess` is used to search for the NASA data collections. Various query parameters can be used to search collections and granules using attributes associated with them in the metadata. More details can be found [here](https://github.com/nsidc/earthaccess/blob/main/notebooks/Demo.ipynb). \n",
    "Below, the CMR Catalog is searched to find collections with `gedi` keywords, managed by `LPCLOUD` provider, and with a `version` number of `002`. The returned response can be used to retrieve the concept-id for each dataset. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96fbc843-1053-4266-bcaf-780f9e22b942",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "collections = earthaccess.collection_query().keyword('gedi').version('002').provider('LPCLOUD').get()\n",
    "pprint(collections[0].summary())\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b6de8f2-68cb-4b7d-87c9-e20f74c0de4e",
   "metadata": {},
   "source": [
    "## Search for GEDI Granules"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a77410fb-ca69-4375-a740-dfca8932f419",
   "metadata": {},
   "source": [
    "Collections `concept-id` is needed to search for data. Below, a dictionary is created to store GEDI V2 collection IDs distributed by LP DAAC, and `GEDI02_B` is selected. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cad25579-3aab-4e9d-ab2b-64dc42b4efc5",
   "metadata": {},
   "outputs": [],
   "source": [
    "gedi_collectionIDs = {}\n",
    "for c in collections:\n",
    "    gedi_collectionIDs[c.summary()['short-name']] = c.summary()['concept-id']\n",
    "gedi_collectionIDs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85e5cf82-1948-427b-975a-472e66ab986f",
   "metadata": {},
   "outputs": [],
   "source": [
    "gedi_products =  ['GEDI02_B']   #['GEDI01_B', 'GEDI02_B', 'GEDI02_A']\n",
    "\n",
    "conceptID = [gedi_collectionIDs[g] for g in gedi_products]\n",
    "conceptID"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e4bd31c-98a7-49d2-905a-dbb60d5f939e",
   "metadata": {},
   "source": [
    "Next, define a temporal range for the query. \n",
    "Please note that GEDI was moved temporarily into hibernation upon the completion of its first mission, which lasted from December 2018 to March 2023. \n",
    "GEDI had been temporarily on pause after March 2023. The instrument returned to the original location on ISS on April 22nd, 2024. GEDI has been collecting data after its return on April 23, 2024, but data has not been publicaly distributed yet. More details can be found [here](https://lpdaac.usgs.gov/news/nasa-announces-pause-in-gedi-mission/).\n",
    "Data is available for the first part of the mission from `2019-04-18` to `2023-03-16` and newly collected data will be available after the proper validations are complete. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03422d51-6052-412b-8884-bda0eb0d504f",
   "metadata": {},
   "outputs": [],
   "source": [
    "tempRange = ('2022-04-01', '2022-05-31') "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7813a2ec-d225-4a40-928e-ffc0af94fce5",
   "metadata": {},
   "source": [
    "A GeoJSON file is used to define the spatial region of interest (ROI). For this example, the ROI is an area in Sequoia National Forest, CA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72c4b8e1-89eb-401e-bf2c-698abaca2a93",
   "metadata": {},
   "outputs": [],
   "source": [
    "polygon = gp.read_file('data/sequoia.geojson')\n",
    "polygon['geometry'][0]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b07d5d6f-899d-41f5-832a-9c3d02d916da",
   "metadata": {},
   "source": [
    "Next, submit the query using `search_data` function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef652862-bf59-4449-b9e3-dad074ac36d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "params = {\n",
    "    \"concept_id\" : conceptID,\n",
    "    \"temporal\": tempRange,\n",
    "    \"polygon\": list(polygon['geometry'][0].exterior.coords),\n",
    "    # bounding_box = bbx,\n",
    "    \"count\": 200\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2912219f-12bf-4b5c-a659-ccda823db7c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "results = earthaccess.search_data(**params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19495596-0aae-42b5-ab12-c0123c664e64",
   "metadata": {},
   "outputs": [],
   "source": [
    "results[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce103519-bfdd-4ac6-b45b-617753371bb2",
   "metadata": {},
   "source": [
    "## Accessing GEDI Data\n",
    "There are two options to access NASA Earth science data stored in the [Earthdata Cloud](https://www.earthdata.nasa.gov/technology/cloud-computing). You can download the data using the `HTTPS` links, create your subset and then work with data locally. The other option is loading data in the memory and only save a subset of data. Loading data in the memory will work both when you are working locally (using `HTTPS` links) or in the cloud (using `S3` links). If you have access to cloud preferably in the same Amazon Web Services (AWS) region us-west2,  you can access and work with data virtually in a cloud-based environment using the `S3` links and skip the downloading part. This method is called “Direct Cloud Access” or, “Direct Access”. Please note that direct access using the `S3` links is only possible if you are working in the Amazon Web Services (AWS) Region us-west-2. If you are working with data locally, you still can load data usig `HTTPS` links but the process could be slower. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4cb861d4-d951-4223-8cd7-03540f0200f2",
   "metadata": {},
   "source": [
    "### Option 1: Donwnloading GEDI data using `HTTPS` links \n",
    "Below, the `HTTPS` links are printed but data can be downloaded using the `download` function from `earthaccess` package directly. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b63ac317-5155-4b24-a356-0a3dd53b6e97",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "data_links = [granule.data_links(access=\"external\") for granule in results]\n",
    "data_links"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66294e21-4b2a-450d-b0f0-aee122b24f7a",
   "metadata": {},
   "source": [
    "Below, the first two granules in the result response is downlaoded below but you can adjust to download them all by replacing `results` with `results[0:2]`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2744cf7f-e42e-43df-b79a-9c972fdcc0f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Only downloaded the first 2 granules\n",
    "# downloaded_files = earthaccess.download(\n",
    "#     results[0:3],\n",
    "#     local_path='data/',\n",
    "# )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8aa15881-2fd9-4e01-8e59-d0f8dec07061",
   "metadata": {},
   "source": [
    "Once your data is downloaded, you can use the **GEDI Subsetter** available in [GEDI-Data_Resources repository](https://github.com/nasa/GEDI-Data-Resources) to subset GEDI sub-orbit granules to your spatial bound and your layers. It is easier to [clone the GEDI-Data-Resources repository](https://github.com/nasa/GEDI-Data-Resources?tab=readme-ov-file#getting-started) and run this notebook but if you have not, you need to adjust the directories below to where the `GEDI_Subsetter.py` is stored. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "557aa0d6-0adc-4d0b-9e9e-070d27beb7e1",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# !python python/scripts/GEDI_Subsetter/GEDI_Subsetter.py --dir data --roi data/sequoia.geojson --beams BEAM0101 --sds '/beam,/quality_flag,/rh,/pai,/pai_z,/pavd_z,/rh100'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f2ce35b-2f19-4491-aa36-12f40a84ec39",
   "metadata": {},
   "source": [
    "As a final step, you can clean up and delete the downloaded source files to free your local space. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38794842-5932-4e27-bd2f-e08735d0fa93",
   "metadata": {},
   "outputs": [],
   "source": [
    "# downloaded = [i for i in os.listdir('data') if i.endswith('.h5')]\n",
    "# for file in downloaded: \n",
    "#     os.remove(f'data/{file}') "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37a84f48-e1c8-4ac5-bd82-f343ee9c024c",
   "metadata": {},
   "source": [
    "### "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a3958f8-283d-45fc-a44a-2b352a467972",
   "metadata": {},
   "source": [
    "### Option 2: accessing GEDI data stored in Earthdata Cloud\n",
    "\n",
    "The other option is loading GEDI data in the memory and only saving a subset of data locally. GEDI V002 data is stored in Earthdata Cloud enabling us to access data in a different way than downloading them. The `open` function from `earthaccess` package will provide a straightforward way to access the data stored in the cloud.  The `open` function offers read and write access for our results to a particular key using a context manager. This will automatically handle the authentication and configurations needed when working locally or in the cloud. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb29ec73-0e25-424c-8569-aca028632e3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "files = earthaccess.open(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be9d2b29",
   "metadata": {},
   "source": [
    "`h5py` package is used to read **GEDI HDF5** GEDI files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd6d58e2-a7ad-4cca-8975-b67ac2288429",
   "metadata": {},
   "outputs": [],
   "source": [
    "gedi_ds = h5py.File(files[0],'r')\n",
    "print(gedi_ds.keys())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bcb43aec-162a-4284-9ed6-f9c8d3d89eb9",
   "metadata": {},
   "outputs": [],
   "source": [
    "gedi_ds['METADATA']['DatasetIdentification'].attrs['shortName']"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10142aca-d2af-41bb-80e2-18acdd113dfe",
   "metadata": {},
   "source": [
    "The available layers (`variables`) and datasets (`datasets`) in a `GEDI_L2B` granule can be accessed (see the commented cell below), but it will take longer for this cell to run. That is why the available GEDI datasets are saved into a JSON file (`GEDI_Datasets.json`) stored in the `data` folder and will be used here. \n",
    "\n",
    "To learn more about the available layers, you can view the GEDI Dictionaries provided in [GEDI products' DOI Landing pages](https://lpdaac.usgs.gov/product_search/?query=gedi&status=Operational&view=cards&sort=title). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b858f798-26d3-40ee-9c5d-8961d97df47a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# variables = []\n",
    "# gedi_ds.visit(variables.append)\n",
    "\n",
    "# datasets = [v.split('/', 1)[-1] for v in variables if isinstance(gedi_ds[v], h5py.Dataset)]\n",
    "# list(set(datasets))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0ff4234",
   "metadata": {},
   "outputs": [],
   "source": [
    "del files, gedi_ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f41e4b59-f63e-44ef-80ac-884488ec6916",
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('data/GEDI_Datasets.json', 'r') as fp:\n",
    "    gedi_var = json.load(fp)\n",
    "\n",
    "gedi_var.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42348796-2fc2-4a45-8852-4ab7f54c1514",
   "metadata": {},
   "source": [
    "Print 'gedi_var' to view the first 20 layers available in `GEDI_L2B` collection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b2d47c1-c004-4166-909a-6b29cf56072b",
   "metadata": {},
   "outputs": [],
   "source": [
    "gedi_var['GEDI_L2B'][0:20]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3cb56d3-4994-4b9c-b8f6-29108d65fe9a",
   "metadata": {},
   "source": [
    "First, Define a subset of datasets you are interested in as alist. If you want to keep all the available datasets save `gedi_var['GEDI_L2B']` in a list. In this example, we defined a subset list for the `GEDI_L2B` product but you can subset more L1 & L2 GEDI products here if your query includes more than one product. To adjust the code, you only need to define subset data in a separate list for each product (`subset_L2A`, `subset_L2B`, and `subset_L1B`) and create a dictionary of all products and subset of layers (see the commented cell below). \n",
    "\n",
    "**Note that the `lat_lowestmode` and `lon_lowestmode` for both GEDI L2A and L2B products are used as reference Latitude and Longitude for each shot. For GEDI L1B, `latitude_bin0` and `longitude_bin0` can be used as reference Latitude and Longitude for each shot. The reference Latitude and Longitude for each shot in addition to `beam` and `shot_number` are included in the subset outputs by default.**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4be15b72-9da7-4565-ae8b-48aaa470d5df",
   "metadata": {},
   "outputs": [],
   "source": [
    "subset_L2B = ['geolocation/degrade_flag', 'geolocation/digital_elevation_model', 'geolocation/elev_lowestmode', 'lat_highestreturn', 'geolocation/lon_highestreturn', 'geolocation/elev_highestreturn', 'l2b_quality_flag', 'rh100', 'pai', 'pai_z', 'pavd_z']\n",
    "\n",
    "subset_data = {'GEDI_L2B': subset_L2B }\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3302fce7-679d-49e1-97ad-f7dd2be5b889",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # example of subsetting data for all three GEDI products distributed by LP DAAC\n",
    "# subset_L2A = gedi_var['GEDI_L2A'] \n",
    "# subset_L2B = gedi_var['GEDI_L2B'] \n",
    "# subset_L1B = gedi_var['GEDI_L1B'] \n",
    "\n",
    "# subset_data = {'GEDI_L2B': subset_L2B, 'GEDI_L2A': subset_L2A, 'GEDI_L1B': subset_L1B }               \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eafa1f38-321e-4428-90ae-99c582d6b4b0",
   "metadata": {},
   "source": [
    "Next, the layers in our subset list provided above is compared with layers stored in `GEDI_Datasets.json`. This step verifies layers are valid and removes the layers that are not. It also creates full path for the layers in the list if that is not provided by user. For example, the dataset 'lat_highestreturn' is in our subset list but the full path should be 'geolocation/lat_highestreturn'. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58383bfb-049d-4f5e-9311-a47c47248b82",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "default = {'GEDI_L1B': ['shot_number', 'beam', 'geolocation/latitude_bin0', 'geolocation/longitude_bin0'],\n",
    "           'GEDI_L2A': ['shot_number', 'beam', 'lat_lowestmode','lon_lowestmode'],\n",
    "           'GEDI_L2B': ['geolocation/shot_number', 'beam', 'geolocation/lat_lowestmode','geolocation/lon_lowestmode']}\n",
    "\n",
    "subset_var = {}\n",
    "\n",
    "for p in list(subset_data.keys()):\n",
    "    subset = []\n",
    "    [subset.append(d) for d in subset_data[p] if d not in subset]\n",
    "    [subset.append(d) for d in default[p] if d not in subset]\n",
    "    datasets_p = []\n",
    "    for s in subset:\n",
    "        my_var = [v for v in gedi_var[p] if v.endswith(f'{s}')]\n",
    "        if len(my_var) == 1:\n",
    "            datasets_p.append(my_var[0])\n",
    "            \n",
    "        elif len(my_var) > 1:\n",
    "            my_var = [v for v in my_var if v.startswith(f'{s}')]\n",
    "                \n",
    "            for l in my_var:\n",
    "                if l not in datasets_p:\n",
    "                    datasets_p.append(l) \n",
    "    \n",
    "    subset_var[p] = datasets_p\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "659b9c5c-5319-4f1e-b94e-18ee73b77e69",
   "metadata": {},
   "outputs": [],
   "source": [
    "subset_var"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce0cafae-317b-43ee-a795-8af193e828fa",
   "metadata": {},
   "source": [
    "In addition to layer subsetting, you can subset layers using specific beams. For instance, you can only select Full Power beams ('BEAM0101', 'BEAM0110', 'BEAM1000', 'BEAM1011')."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d36bb79-3e22-47e3-80b9-c78a80c871bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "beams = ['BEAM0000', 'BEAM0001', 'BEAM0010', 'BEAM0011', 'BEAM0101', 'BEAM0110', 'BEAM1000', 'BEAM1011'] # ['BEAM0101', 'BEAM0110', 'BEAM1000', 'BEAM1011']\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd1b3a27-ce27-4b78-b69a-34cc5c0a210a",
   "metadata": {},
   "source": [
    "Below, functions are defined to subset the GEDI granules using the beams, datasets, and area of interest. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf23f6bb-237c-4a13-b1aa-bde30ee17e52",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gedi_to_dataframe(granule, beams, vars):\n",
    "    \"\"\"\n",
    "    This function takes existing method of getting data from a GEDI hdf5\n",
    "    and makes it dynamic so it will retrieve subset of beams and variables from a list provided by the user.\n",
    "    All column names are taken from the hdf5 source file.\n",
    "    \"\"\"\n",
    "\n",
    "    ds = earthaccess.open([granule])[0]\n",
    "    #read the dataset\n",
    "    gedi_ds = h5py.File(ds,'r')\n",
    "    # see what is the data product \n",
    "    product = gedi_ds['METADATA']['DatasetIdentification'].attrs['shortName']\n",
    "    print(product)\n",
    "    fileName = gedi_ds['METADATA']['DatasetIdentification'].attrs['fileName']\n",
    "    date = datetime.strptime(fileName.rsplit('_')[2], '%Y%j%H%M%S').strftime('%Y-%m-%d %H:%M:%S')\n",
    "    # Create empty DataFrame for this beam\n",
    "    df_beam = pd.DataFrame(columns=vars[product])\n",
    "    \n",
    "    for b in beams:\n",
    "        data_dic = {}\n",
    "        for v in vars[product]:\n",
    "            # print(b,v)\n",
    "            value = gedi_ds[f'{b}/{v}'][()]\n",
    "            data_dic[v] = value.tolist() \n",
    "            \n",
    "        df_beam = pd.concat([df_beam, pd.DataFrame(data_dic)],join=\"inner\")\n",
    "        \n",
    "        # add product, beam, file name, and date columns \n",
    "        df_beam.insert(0, 'product', product)\n",
    "        df_beam.insert(1, 'Beam', b)\n",
    "        df_beam.insert(2, 'fileName' , fileName)\n",
    "        df_beam.insert(3, 'date', date)\n",
    "\n",
    "    # rename the latitude and longitude here to simplify the dataframe\n",
    "    df_beam= df_beam.rename(columns={'geolocation/lat_lowestmode': 'lat', 'geolocation/lon_lowestmode': 'lon', \n",
    "                                     'geolocation/latitude_bin0':'lat_bin0', 'geolocation/longitude_bin0':'lon_bin0',\n",
    "                                     'geolocation/shot_number':'shot_number',\n",
    "                                     'lat_lowestmode': 'lat', 'lon_lowestmode': 'lon'})\n",
    "\n",
    "    return(product, df_beam.reset_index(drop=True))\n",
    "\n",
    "\n",
    "def clip_gedi(dataframe,geojson):\n",
    "    \"\"\"\n",
    "    This function takes the subset of GEDI data stored in a Geopandas dataframe and creates a spatial subset.\n",
    "    \"\"\"\n",
    "    #read the GeoJSON\n",
    "    ROI = gp.GeoDataFrame.from_file(geojson)\n",
    "    ROI.crs = 'EPSG:4326'\n",
    "\n",
    "    # Take the lat/lon dataframe and convert each lat/lon to a shapely point and convert to a Geodataframe\n",
    "    try:\n",
    "        dataframe = gp.GeoDataFrame(dataframe, geometry=dataframe.apply(lambda row: Point(row.lon, row.lat), axis=1))\n",
    "    except:\n",
    "        dataframe = gp.GeoDataFrame(dataframe, geometry=dataframe.apply(lambda row: Point(row.lon_bin0, row.lat_bin0), axis=1))\n",
    "        \n",
    "    dataframe = dataframe.set_crs('EPSG:4326')\n",
    "\n",
    "    shot_list = []\n",
    "    for num, geom in enumerate(dataframe['geometry']):\n",
    "        if ROI.contains(geom)[0]:\n",
    "            shot_n = dataframe.loc[num, 'shot_number']\n",
    "            shot_list.append(shot_n)\n",
    "            \n",
    "\n",
    "    DF = dataframe.where(dataframe['shot_number'].isin(shot_list))\n",
    "    DF = DF.dropna().reset_index(drop=True)\n",
    "    return DF"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adc59f87-7b4a-4eb6-94e9-fe73c15933d4",
   "metadata": {},
   "source": [
    "A separate empty dataframe is created for each product. In the first step, source granules are accessed and a subset of data by layer is created.  Next, the spatial subset is implemented and the output Geodataframe object is returned. Finally, the subsets for all the source files are concatenated into a single Geodataframe. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "c53d5f8f-4039-4aff-b575-87d3d1d679f7",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "l2a_df = pd.DataFrame()\n",
    "l2b_df = pd.DataFrame()\n",
    "l1b_df = pd.DataFrame()\n",
    "num = 0\n",
    "for granule in results:\n",
    "    num += 0\n",
    "    # print(granule)\n",
    "    # subset data by band\n",
    "    prod, subset_df = gedi_to_dataframe(granule,beams,subset_var)\n",
    "    # clip to ROI    \n",
    "    subset_df_clip = clip_gedi(subset_df, 'data/sequoia.geojson')\n",
    "    if 'L2A' in prod:\n",
    "        l2a_df = gp.GeoDataFrame(pd.concat([l2a_df, subset_df_clip]))\n",
    "    elif 'L2B' in prod:\n",
    "        l2b_df = gp.GeoDataFrame(pd.concat([l2b_df, subset_df_clip]))\n",
    "    elif 'L1B' in prod:\n",
    "        l1b_df = gp.GeoDataFrame(pd.concat([l1b_df, subset_df_clip]))\n",
    "\n",
    "    del subset_df_clip, subset_df, prod\n",
    "\n",
    "# Reset the indeces\n",
    "l1b_df = l1b_df.reset_index(drop=True)\n",
    "l2a_df = l2a_df.reset_index(drop=True)\n",
    "l2b_df = l2b_df.reset_index(drop=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cab89006-559d-4cea-a2ab-fa40a932d061",
   "metadata": {},
   "outputs": [],
   "source": [
    "l1b_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d42f5ae-773a-4ede-831a-31c6c62657b0",
   "metadata": {},
   "source": [
    "Data types are adjusted to save the GeoJSON file successfully. In this example, `pai_z` and `pavd_z` data types are updated to string. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2764a76-ff25-474e-af69-b882409d8973",
   "metadata": {},
   "outputs": [],
   "source": [
    "l2b_df.columns\n",
    "for c in ['pai_z', 'pavd_z']:\n",
    "    l2b_df[c] =  str(l2b_df[c])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "897e537a-884e-4064-ab44-b26d5b789770",
   "metadata": {},
   "source": [
    "Finally, we exported subset files as a GeoJSON using the `GeoDataFrame.to_file` function from `geopandas` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "185daa55-28ab-41ef-b905-6cb9350436f4",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from datetime import date\n",
    "today = date.today()\n",
    "\n",
    "if len(l1b_df) != 0:\n",
    "    l1b_df.to_file(f'data/GEDI_L1B_{today}.geojson', driver='GeoJSON') \n",
    "\n",
    "if len(l2a_df) != 0:\n",
    "    l2a_df.to_file(f'data/GEDI_L2A_{today}.geojson', driver='GeoJSON') \n",
    "\n",
    "if len(l2b_df) != 0:\n",
    "    l2b_df.to_file(f'data/GEDI_L2B_{today}.geojson', driver='GeoJSON')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f0d679e-1648-4724-b5b6-ab1dd1d64ca2",
   "metadata": {},
   "source": [
    "The GeoJSN object can be viewed using QGIS or programmatically. Below, a simple map is created to visualize the spatial coverage of data over our ROI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2a40eaf-a94a-4602-b039-de4bb7aa81d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from geoviews import opts, tile_sources as gvts\n",
    "import geoviews\n",
    "geoviews.extension('bokeh')\n",
    "\n",
    "gvts.EsriImagery * geoviews.Points(l2b_df, vdims=['date']).options(tools=['hover'], \n",
    "                                                                   height=900, width=900, size=1, color='yellow', \n",
    "                                                                   fontsize={'xticks': 10, 'yticks': 10, 'xlabel':16, 'ylabel': 16})\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b751275c-b79e-48cc-a006-4d7ac8a26ec3",
   "metadata": {},
   "outputs": [],
   "source": [
    "l2b_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06e514c2-89de-40c6-a704-4035a4b5d18c",
   "metadata": {},
   "source": [
    "\n",
    "## Contact Info:  \n",
    "\n",
    "Email: LPDAAC@usgs.gov  \n",
    "Voice: +1-866-573-3222  \n",
    "Organization: Land Processes Distributed Active Archive Center (LP DAAC)¹  \n",
    "Website: <https://lpdaac.usgs.gov/>  \n",
    "Date last modified: 02-20-2024  \n",
    "\n",
    "¹Work performed under USGS contract G15PD00467 for NASA contract NNG14HH33I.  "
   ]
  }
 ],
 "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.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}