notebooks/Download_GEDI_L4A.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d02c930f",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-06-08T16:09:52.682153Z",
     "iopub.status.busy": "2024-06-08T16:09:52.681797Z",
     "iopub.status.idle": "2024-06-08T16:14:34.501167Z",
     "shell.execute_reply": "2024-06-08T16:14:34.499948Z"
    },
    "id": "pDJ2OijZX5Kk",
    "outputId": "85b9fd06-140d-437f-925f-d0c52eef8692",
    "papermill": {
     "duration": 281.828502,
     "end_time": "2024-06-08T16:14:34.503739",
     "exception": false,
     "start_time": "2024-06-08T16:09:52.675237",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "!pip install contextily~=1.3.0\n",
    "!pip install geopandas~=0.13.2\n",
    "!pip install folium~=0.14.0\n",
    "!pip install h5py~=3.8.0\n",
    "!pip install harmony-py~=0.4.9\n",
    "!pip install IPython~=8.12.2\n",
    "!pip install matplotlib~=3.7.1\n",
    "!pip install matplotlib_scalebar~=0.8.1\n",
    "!pip install netCDF4~=1.6.4\n",
    "!pip install numpy~=1.24.3\n",
    "!pip install owslib~=0.29.2\n",
    "!pip install pandas~=2.0.2\n",
    "!pip install pydap~=3.4.1\n",
    "!pip install rasterio~=1.3.7\n",
    "!pip install rasterstats~=0.19.0\n",
    "!pip install requests~=2.31.0\n",
    "!pip install s3fs~=2023.6.0\n",
    "!pip install seaborn~=0.12.2\n",
    "!pip install scipy~=1.10.1\n",
    "!pip install shapely~=2.0.1\n",
    "!pip install tabulate~=0.9.0\n",
    "!pip install mlxtend\n",
    "!pip install plotly\n",
    "!pip install scikit-learn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "38881f92",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-06-08T16:14:34.539257Z",
     "iopub.status.busy": "2024-06-08T16:14:34.538310Z",
     "iopub.status.idle": "2024-06-08T16:14:36.960672Z",
     "shell.execute_reply": "2024-06-08T16:14:36.959766Z"
    },
    "id": "SgjKuqsuYr1e",
    "papermill": {
     "duration": 2.442547,
     "end_time": "2024-06-08T16:14:36.963193",
     "exception": false,
     "start_time": "2024-06-08T16:14:34.520646",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# import modules\n",
    "%matplotlib inline\n",
    "import os\n",
    "import traceback\n",
    "from os import path\n",
    "from glob import glob\n",
    "import requests\n",
    "import datetime as dt\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import contextily as ctx\n",
    "import netCDF4 as nc\n",
    "from matplotlib import pyplot as plt\n",
    "from pydap.cas.urs import setup_session\n",
    "from shapely.ops import orient\n",
    "from IPython import display\n",
    "from requests.adapters import HTTPAdapter, Retry\n",
    "import seaborn as sns\n",
    "sns.set(style='whitegrid')\n",
    "pd.set_option('mode.chained_assignment', None)\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67242884",
   "metadata": {},
   "source": [
    "#### Add the folder path of where the .shp resides in `input_shape_file_path` variable\n",
    "#### Add the filename of .shp file in `input_shape_file_name` variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b1e7d9b1",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-06-08T16:14:36.998164Z",
     "iopub.status.busy": "2024-06-08T16:14:36.997628Z",
     "iopub.status.idle": "2024-06-08T16:14:37.003720Z",
     "shell.execute_reply": "2024-06-08T16:14:37.002873Z"
    },
    "id": "JT17HQJ05LR5",
    "papermill": {
     "duration": 0.025966,
     "end_time": "2024-06-08T16:14:37.005820",
     "exception": false,
     "start_time": "2024-06-08T16:14:36.979854",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "input_shape_file_path = \"/kaggle/input/nc-forests-all-time\"\n",
    "input_shape_file_name = \"NC_forests_all_time.shp\"\n",
    "os.makedirs(input_shape_file_path, exist_ok=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80c3e4ee",
   "metadata": {},
   "source": [
    "#### Add the year you are interested in `year_` variable (We downloaded data between 2018 - 2021)\n",
    "#### Add the folder path where **Gedi L4A** should be downloaded in `outdir_folder` variable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ffa21442",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-06-08T16:14:37.040252Z",
     "iopub.status.busy": "2024-06-08T16:14:37.039594Z",
     "iopub.status.idle": "2024-06-08T16:14:37.044912Z",
     "shell.execute_reply": "2024-06-08T16:14:37.043988Z"
    },
    "id": "8uKGEay5YsFa",
    "papermill": {
     "duration": 0.02471,
     "end_time": "2024-06-08T16:14:37.046951",
     "exception": false,
     "start_time": "2024-06-08T16:14:37.022241",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Creating the path to the GEDI data\n",
    "year_ = 2021\n",
    "outdir_folder = \"/kaggle/working/BulkDownload\"\n",
    "outdir_output = os.path.join(outdir_folder, \"agbd_output/{}\".format(year_))\n",
    "os.makedirs(outdir_output, exist_ok=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "06a790f8",
   "metadata": {},
   "source": [
    "#### Add your Earthdata **username** and **password**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ceda21c2",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-06-08T16:14:37.082958Z",
     "iopub.status.busy": "2024-06-08T16:14:37.082238Z",
     "iopub.status.idle": "2024-06-08T16:14:37.086442Z",
     "shell.execute_reply": "2024-06-08T16:14:37.085438Z"
    },
    "id": "gQ-_chm7Yw5T",
    "papermill": {
     "duration": 0.024935,
     "end_time": "2024-06-08T16:14:37.088348",
     "exception": false,
     "start_time": "2024-06-08T16:14:37.063413",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Setting up the Earthdata URS Login session.\n",
    "# Get account here https://urs.earthdata.nasa.gov/\n",
    "# Replace username and password with your credentials\n",
    "\n",
    "username = \"MyEarthDataUsername\"\n",
    "password= \"MyEarthDataPassword\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "63016a08",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-06-08T16:14:37.122779Z",
     "iopub.status.busy": "2024-06-08T16:14:37.122371Z",
     "iopub.status.idle": "2024-06-08T16:14:37.127221Z",
     "shell.execute_reply": "2024-06-08T16:14:37.126432Z"
    },
    "id": "mhX-6w1tY-Lr",
    "papermill": {
     "duration": 0.024495,
     "end_time": "2024-06-08T16:14:37.129247",
     "exception": false,
     "start_time": "2024-06-08T16:14:37.104752",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# variables of interests\n",
    "variables = ['agbd', 'l4_quality_flag', 'degrade_flag', 'agbd_se', 'land_cover_data/pft_class', 'solar_elevation']\n",
    "\n",
    "# beams of interest\n",
    "beams = ['BEAM0000', 'BEAM0001', 'BEAM0010', 'BEAM0011', 'BEAM0101', 'BEAM0110', 'BEAM1000', 'BEAM1011']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "7ced7307",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-06-08T16:14:37.163635Z",
     "iopub.status.busy": "2024-06-08T16:14:37.163143Z",
     "iopub.status.idle": "2024-06-08T16:14:37.167688Z",
     "shell.execute_reply": "2024-06-08T16:14:37.166741Z"
    },
    "id": "kzBF7FqsZUhM",
    "papermill": {
     "duration": 0.024232,
     "end_time": "2024-06-08T16:14:37.169804",
     "exception": false,
     "start_time": "2024-06-08T16:14:37.145572",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "headers = ['lat_lowestmode', 'lon_lowestmode', 'elev_lowestmode', 'shot_number']\n",
    "headers.extend(variables)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "0b09db1c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-06-08T16:14:37.238080Z",
     "iopub.status.busy": "2024-06-08T16:14:37.237300Z",
     "iopub.status.idle": "2024-06-08T16:14:37.246926Z",
     "shell.execute_reply": "2024-06-08T16:14:37.245988Z"
    },
    "id": "u3CF9jnAY_DV",
    "papermill": {
     "duration": 0.029674,
     "end_time": "2024-06-08T16:14:37.248968",
     "exception": false,
     "start_time": "2024-06-08T16:14:37.219294",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def search_all_data(aca, index_):\n",
    "  doi = '10.3334/ORNLDAAC/2056'\n",
    "  # CMR API base url\n",
    "  cmrurl='https://cmr.earthdata.nasa.gov/search/'\n",
    "  doisearch = f\"{cmrurl}collections.json?doi={doi}\"\n",
    "  concept_id = requests.get(doisearch).json()['feed']['entry'][0]['id']\n",
    "\n",
    "  # defining geojson\n",
    "\n",
    "\n",
    "  geojson = {\"shapefile\": (\"aca.geojson\", aca.geometry.to_json(), \"application/geo+json\")}\n",
    "  page_num = 1\n",
    "  page_size = 2000 # CMR page size limit. do not change this. it is maxed out\n",
    "\n",
    "  start_date = dt.datetime(year_, 1, 1) # specify your own start date\n",
    "  end_date = dt.datetime(year_, 12, 31)  # specify your end start date\n",
    "\n",
    "  # CMR formatted start and end times\n",
    "  dt_format = '%Y-%m-%dT%H:%M:%SZ'\n",
    "  temporal_str = start_date.strftime(dt_format) + ',' + end_date.strftime(dt_format)\n",
    "\n",
    "  opendap_arr = []\n",
    "\n",
    "  while True:\n",
    "\n",
    "      # defining parameters\n",
    "      cmr_param = {\n",
    "          \"collection_concept_id\": concept_id,\n",
    "          \"temporal\": temporal_str,\n",
    "          \"page_size\": page_size,\n",
    "          \"page_num\": page_num,\n",
    "          \"simplify-shapefile\": 'true' # this is needed to bypass 5000 coordinates limit of CMR\n",
    "      }\n",
    "\n",
    "      granulesearch = f\"{cmrurl}granules.json\"\n",
    "      response = requests.post(granulesearch, data=cmr_param, files=geojson)\n",
    "      granules = response.json()['feed']['entry']\n",
    "      if granules:\n",
    "          for g in granules:\n",
    "              # Get OPeNDAP URLs\n",
    "              for links in g['links']:\n",
    "                  if 'title' in links and links['title'].startswith('OPeNDAP'):\n",
    "                      opendap_url = links['href']\n",
    "                      opendap_arr.append(opendap_url)\n",
    "          page_num += 1\n",
    "      else:\n",
    "          break\n",
    "\n",
    "  total_granules = len(opendap_arr)\n",
    "  print (\"Total granules found: \", total_granules)\n",
    "\n",
    "  return {index_: opendap_arr}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "0fa4aa5a",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-06-08T16:14:37.283749Z",
     "iopub.status.busy": "2024-06-08T16:14:37.283010Z",
     "iopub.status.idle": "2024-06-08T16:14:37.296198Z",
     "shell.execute_reply": "2024-06-08T16:14:37.295255Z"
    },
    "id": "h0adJCgRZWh7",
    "papermill": {
     "duration": 0.032767,
     "end_time": "2024-06-08T16:14:37.298188",
     "exception": false,
     "start_time": "2024-06-08T16:14:37.265421",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def download_data(opendap_map, aca):\n",
    "  session = setup_session(username, password, check_url=\"https://opendap.earthdata.nasa.gov/\")\n",
    "  index_ = list(opendap_map.keys())[0]\n",
    "  opendap_arr = opendap_map[index_]\n",
    "  csv_file_name = \"{}.csv\".format(index_)\n",
    "  out_csv = os.path.join(outdir_output,csv_file_name)\n",
    "\n",
    "  with open(out_csv, \"w\") as f:\n",
    "      f.write(','.join(headers)+'\\n')\n",
    "  retries = Retry(total=50, backoff_factor=0.1, status_forcelist=[ 500, 502, 503, 504 ])\n",
    "  session.mount('https://', HTTPAdapter(max_retries=retries))\n",
    "\n",
    "  c=0\n",
    "  for g_name in opendap_arr:\n",
    "      c += 1\n",
    "      # loop over all beams\n",
    "      for beam in beams:\n",
    "          # 1. Retrieving lat, lon coordinates for the file\n",
    "          hyrax_url = f\"{g_name}.dap.nc4?dap4.ce=/{beam}/lon_lowestmode;/{beam}/lat_lowestmode\"\n",
    "          r = session.get(hyrax_url)\n",
    "          if (r.status_code != 400):\n",
    "              ds = nc.Dataset('hyrax', memory=r.content)\n",
    "              lat = ds[beam]['lat_lowestmode'][:]\n",
    "              lon = ds[beam]['lon_lowestmode'][:]\n",
    "              ds.close()\n",
    "              df = pd.DataFrame({'lat_lowestmode': lat, 'lon_lowestmode': lon}) # creating pandas dataframe\n",
    "\n",
    "              # 2. Subsetting by bounds of the area of interest\n",
    "              # converting to geopandas dataframe\n",
    "              gdf = gpd.GeoDataFrame(df, geometry=gpd.points_from_xy(df.lon_lowestmode, df.lat_lowestmode))\n",
    "              gdf_aca = gdf[gdf['geometry'].within(aca.geometry[index_])]\n",
    "              if not gdf_aca.empty:\n",
    "                  # creating empty columns for variables\n",
    "                  for v in headers[2:]:\n",
    "                      gdf_aca[v] = None\n",
    "                  # 3. retrieving variables of interest, agbd, agbd_t in this case.\n",
    "                  # We are only retriving the shots within subset area.\n",
    "                  for _, df_gr in gdf_aca.groupby((gdf_aca.index.to_series().diff() > 1).cumsum()):\n",
    "                      i = df_gr.index.min()\n",
    "                      j = df_gr.index.max()\n",
    "                      for v in headers[2:]:\n",
    "                          var_s = f\"/{beam}/{v}%5B{i}:{j}%5D\"\n",
    "                          hyrax_url = f\"{g_name}.dap.nc4?dap4.ce={var_s}\"\n",
    "                          r = session.get(hyrax_url)\n",
    "                          if (r.status_code != 400):\n",
    "                              ds = nc.Dataset('hyrax', memory=r.content)\n",
    "                              gdf_aca.loc[i:j, (v)] = ds[beam][v][:]\n",
    "                              ds.close()\n",
    "\n",
    "                  # saving the output file\n",
    "                  gdf_aca.to_csv(out_csv, mode='a', index=False, header=False, columns=headers)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "eb83d6ed",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-06-08T16:14:37.332861Z",
     "iopub.status.busy": "2024-06-08T16:14:37.332012Z",
     "iopub.status.idle": "2024-06-08T16:14:37.338241Z",
     "shell.execute_reply": "2024-06-08T16:14:37.337316Z"
    },
    "id": "KGFFT7nCZrO7",
    "papermill": {
     "duration": 0.025552,
     "end_time": "2024-06-08T16:14:37.340135",
     "exception": false,
     "start_time": "2024-06-08T16:14:37.314583",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def search_and_download(shape_file_data_frame):\n",
    "  try:\n",
    "    df = shape_file_data_frame[0]\n",
    "    index_ = shape_file_data_frame[1]\n",
    "    new_df = df.copy()\n",
    "    new_df.crs = \"EPSG:4326\"\n",
    "    new_df.geometry = new_df.geometry.apply(orient, args=(1,))\n",
    "    opendap_map = search_all_data(new_df, index_)\n",
    "    download_data(opendap_map, new_df)\n",
    "    print(\"Finished downloading {}\".format(index_))\n",
    "    return 0\n",
    "  except Exception:\n",
    "    # there will definitely be errors while downloading from opendap servers. Ignore them\n",
    "    print(\"Error in downloading {}\".format(index_))\n",
    "    traceback.print_exc()\n",
    "    return 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "e3c55bbc",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-06-08T16:14:37.374777Z",
     "iopub.status.busy": "2024-06-08T16:14:37.373899Z",
     "iopub.status.idle": "2024-06-08T16:14:42.677068Z",
     "shell.execute_reply": "2024-06-08T16:14:42.676190Z"
    },
    "id": "-PhX_3pnZAHY",
    "papermill": {
     "duration": 5.323097,
     "end_time": "2024-06-08T16:14:42.679589",
     "exception": false,
     "start_time": "2024-06-08T16:14:37.356492",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# read the data from .shp file. Essentialy we want to download L4A data for the regions in the shape file\n",
    "all_data = gpd.read_file(os.path.join(input_shape_file_path, input_shape_file_name))\n",
    "all_data.crs = \"EPSG:4326\"\n",
    "all_data.geometry = all_data.geometry.apply(orient, args=(1,))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55146dc0",
   "metadata": {},
   "source": [
    "#### Define the number of processes to use in order to download data from opendap in `no_of_prc` variable\n",
    "#### Every patch/region from the shape file will be a seperate process/download request to the opendap server\n",
    "#### **Note - This download will take a while (days) since we need to download data for every patch. Hence (what we did) we split this task and created multiple notebooks/user accounts which download different patches by manipulating `start` and `end` variable**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f54e7a7e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-06-08T16:14:42.715891Z",
     "iopub.status.busy": "2024-06-08T16:14:42.715512Z",
     "iopub.status.idle": "2024-06-08T16:14:45.973946Z",
     "shell.execute_reply": "2024-06-08T16:14:45.972864Z"
    },
    "id": "1dPLa0R5aG3y",
    "papermill": {
     "duration": 3.279198,
     "end_time": "2024-06-08T16:14:45.976121",
     "exception": false,
     "start_time": "2024-06-08T16:14:42.696923",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "from multiprocessing import Pool\n",
    "no_of_prc = 20\n",
    "start = 0\n",
    "end = len(all_data)\n",
    "matrix = [(all_data.iloc[i:i+1].copy(), i) for i in range(start, end)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a47d61e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-06-08T16:14:46.010911Z",
     "iopub.status.busy": "2024-06-08T16:14:46.010534Z",
     "iopub.status.idle": "2024-06-09T01:12:56.456118Z",
     "shell.execute_reply": "2024-06-09T01:12:56.454360Z"
    },
    "id": "pgKxipuEdeR2",
    "outputId": "85969003-59fa-406f-ccbf-af265923a468",
    "papermill": {
     "duration": 32290.535975,
     "end_time": "2024-06-09T01:12:56.528762",
     "exception": false,
     "start_time": "2024-06-08T16:14:45.992787",
     "status": "completed"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "if __name__ == '__main__':\n",
    "  with Pool(no_of_prc) as p:\n",
    "      p.map(search_and_download, matrix)"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kaggle": {
   "accelerator": "none",
   "dataSources": [
    {
     "datasetId": 5098016,
     "sourceId": 8535108,
     "sourceType": "datasetVersion"
    }
   ],
   "dockerImageVersionId": 30698,
   "isGpuEnabled": false,
   "isInternetEnabled": true,
   "language": "python",
   "sourceType": "notebook"
  },
  "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.10.13"
  },
  "papermill": {
   "default_parameters": {},
   "duration": 32589.748405,
   "end_time": "2024-06-09T01:12:59.887968",
   "environment_variables": {},
   "exception": null,
   "input_path": "__notebook__.ipynb",
   "output_path": "__notebook__.ipynb",
   "parameters": {},
   "start_time": "2024-06-08T16:09:50.139563",
   "version": "2.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}