python/how-tos/How_to_find_EMIT_data_using_CMR_API.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to: Find and Access EMIT Data\n",
    "\n",
    "**Summary**  \n",
    "\n",
    "This notebook will explain how to access Earth Surface Minteral Dust Source Investigation (EMIT) data programmaticly using NASA's CMR API. The Common Metadata Repository (CMR) is a metadata system that catalogs Earth Science data and associated metadata records. The CMR Application Programming Interface (API) provides programatic search capabilities through CMR's vast metadata holdings using various parameters and keywords. When querying NASA's CMR, there is a limit of 1 million granules matched and only 2000 granules returned per page. \n",
    "\n",
    "**Requirements:**\n",
    "+ A NASA [Earthdata Login](https://urs.earthdata.nasa.gov/) account is required to download EMIT data   \n",
    "+ Selected the `emit_tutorials` environment as the kernel for this notebook.\n",
    "  + For instructions on setting up the environment, follow the the `setup_instructions.md` included in the `/setup/` folder of the repository.  \n",
    "\n",
    "**Learning Objectives**  \n",
    "- How to find EMIT data using NASA's CMR API\n",
    "- How to download programmatically "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the required packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import requests\n",
    "import earthaccess\n",
    "import pandas as pd\n",
    "import datetime as dt\n",
    "import geopandas\n",
    "from shapely.geometry import MultiPolygon, Polygon, box"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Obtaining the Concept ID\n",
    "\n",
    "NASA EarthData's unique ID for this dataset (called Concept ID) is needed for searching the dataset. The dataset Digital Object Identifier or DOI can be used to obtain the Concept ID. DOIs can be found by clicking the `Citation` link on the LP DAAC's [EMIT Product Pages](https://lpdaac.usgs.gov/product_search/?query=emit&view=cards&sort=title)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "doi = '10.5067/EMIT/EMITL2ARFL.001'# EMIT L2A Reflectance\n",
    "\n",
    "# CMR API base url\n",
    "cmrurl='https://cmr.earthdata.nasa.gov/search/' \n",
    "\n",
    "doisearch = cmrurl + 'collections.json?doi=' + doi\n",
    "concept_id = requests.get(doisearch).json()['feed']['entry'][0]['id']\n",
    "print(concept_id)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the unique NASA-given concept ID for the EMIT L2A Reflectance dataset, which can be used to retrieve relevant files (or granules)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Searching using CMR API\n",
    "\n",
    "When searching the CMR API, users can provide spatial bounds and date-time ranges to narrow their search. These spatial bounds can be either, points, a bounding box, or a polygon. \n",
    "\n",
    "Specify start time and dates and reformat them to the structure necessary for searching CMR."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Temporal Bound - Year, month, day. Hour, minutes, and seconds (ZULU) can also be included \n",
    "start_date = dt.datetime(2022, 9, 3)\n",
    "end_date = dt.datetime(2022, 9, 3, 23, 23, 59)  \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",
    "print(temporal_str)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The CMR API only allows 2000 results to be shown at a time. Using `page_num` allows a user to loop through the search result pages. The sections below walk through using Points, Bounding Boxes, and Polygons to spatially constrain a search made using the CMR API. \n",
    "\n",
    "### Search using Points\n",
    "\n",
    "To search using a point we specify a latitude and longitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Search using a Point\n",
    "\n",
    "lon = -62.1123\n",
    "lat = -39.89402\n",
    "point_str = str(lon) +','+ str(lat)\n",
    "\n",
    "page_num = 1\n",
    "page_size = 2000 # CMR page size limit\n",
    "\n",
    "granule_arr = []\n",
    "\n",
    "while True:\n",
    "    \n",
    "     # defining parameters\n",
    "    cmr_param = {\n",
    "        \"collection_concept_id\": concept_id, \n",
    "        \"page_size\": page_size,\n",
    "        \"page_num\": page_num,\n",
    "        \"temporal\": temporal_str,\n",
    "        \"point\":point_str\n",
    "    }\n",
    "\n",
    "    granulesearch = cmrurl + 'granules.json'\n",
    "    response = requests.post(granulesearch, data=cmr_param)\n",
    "    granules = response.json()['feed']['entry']\n",
    "       \n",
    "    if granules:\n",
    "        for g in granules:\n",
    "            granule_urls = ''\n",
    "            granule_poly = ''\n",
    "                       \n",
    "            # read cloud cover\n",
    "            cloud_cover = g['cloud_cover']\n",
    "    \n",
    "            # reading bounding geometries\n",
    "            if 'polygons' in g:\n",
    "                polygons= g['polygons']\n",
    "                multipolygons = []\n",
    "                for poly in polygons:\n",
    "                    i=iter(poly[0].split (\" \"))\n",
    "                    ltln = list(map(\" \".join,zip(i,i)))\n",
    "                    multipolygons.append(Polygon([[float(p.split(\" \")[1]), float(p.split(\" \")[0])] for p in ltln]))\n",
    "                granule_poly = MultiPolygon(multipolygons)\n",
    "            \n",
    "            # Get https URLs to .nc files and exclude .dmrpp files\n",
    "            granule_urls = [x['href'] for x in g['links'] if 'https' in x['href'] and '.nc' in x['href'] and '.dmrpp' not in x['href']]\n",
    "            # Add to list\n",
    "            granule_arr.append([granule_urls, cloud_cover, granule_poly])\n",
    "                           \n",
    "        page_num += 1\n",
    "    else: \n",
    "        break\n",
    "print(granule_arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Search using a bounding box\n",
    "For this we'll use a bounding box along the coast of Argentina with a bottom left corner of -62.1123 Longitude, -39.89402 Latitude, and a top right corner of -61.70801 Longitude and -39.57769 Latitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Search Using a Bounding Box\n",
    "bound = (-62.1123, -39.89402, -61.70801, -39.57769) \n",
    "bound_str = ','.join(map(str,bound))\n",
    "\n",
    "page_num = 1\n",
    "page_size = 2000 # CMR page size limit\n",
    "\n",
    "granule_arr = []\n",
    "\n",
    "while True:\n",
    "    \n",
    "     # defining parameters\n",
    "    cmr_param = {\n",
    "        \"collection_concept_id\": concept_id, \n",
    "        \"page_size\": page_size,\n",
    "        \"page_num\": page_num,\n",
    "        \"temporal\": temporal_str,\n",
    "        \"bounding_box[]\":bound_str\n",
    "    }\n",
    "\n",
    "    granulesearch = cmrurl + 'granules.json'\n",
    "    response = requests.post(granulesearch, data=cmr_param)\n",
    "    granules = response.json()['feed']['entry']\n",
    "       \n",
    "    if granules:\n",
    "        for g in granules:\n",
    "            granule_urls = ''\n",
    "            granule_poly = ''\n",
    "                       \n",
    "            # read cloud cover\n",
    "            cloud_cover = g['cloud_cover']\n",
    "    \n",
    "            # reading results bounding geometries\n",
    "            if 'polygons' in g:\n",
    "                polygons= g['polygons']\n",
    "                multipolygons = []\n",
    "                for poly in polygons:\n",
    "                    i=iter(poly[0].split (\" \"))\n",
    "                    ltln = list(map(\" \".join,zip(i,i)))\n",
    "                    multipolygons.append(Polygon([[float(p.split(\" \")[1]), float(p.split(\" \")[0])] for p in ltln]))\n",
    "                granule_poly = MultiPolygon(multipolygons)\n",
    "            \n",
    "            # Get https URLs to .nc files and exclude .dmrpp files\n",
    "            granule_urls = [x['href'] for x in g['links'] if 'https' in x['href'] and '.nc' in x['href'] and '.dmrpp' not in x['href']]\n",
    "            # Add to list\n",
    "            granule_arr.append([granule_urls, cloud_cover, granule_poly])\n",
    "                           \n",
    "        page_num += 1\n",
    "    else: \n",
    "        break\n",
    "print(granule_arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Search a Polygon\n",
    "\n",
    "A polygon can also be used to spatially search using the CMR API. A shapefile, geojson, or other format can be opened as a geopandas dataframe, then reformatted to a geojson format to be sent as a parameter in the CMR search. Note that very complex shapefiles must be simplified, there is a 5000 coordinate limit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Search using a Polygon\n",
    "polygon = geopandas.read_file('../../data/isla_gaviota.geojson')\n",
    "geojson = {\"shapefile\": (\"isla_gaviota.geojson\", polygon.geometry.to_json(), \"application/geo+json\")}\n",
    "\n",
    "page_num = 1\n",
    "page_size = 2000 # CMR page size limit\n",
    "\n",
    "granule_arr = []\n",
    "\n",
    "while True:\n",
    "    \n",
    "     # defining parameters\n",
    "    cmr_param = {\n",
    "        \"collection_concept_id\": concept_id, \n",
    "        \"page_size\": page_size,\n",
    "        \"page_num\": page_num,\n",
    "        \"temporal\": temporal_str,\n",
    "        \"simplify-shapefile\": 'true' # this is needed to bypass 5000 coordinates limit of CMR\n",
    "    }\n",
    "\n",
    "    granulesearch = cmrurl + 'granules.json'\n",
    "    response = requests.post(granulesearch, data=cmr_param, files=geojson)\n",
    "    granules = response.json()['feed']['entry']\n",
    "       \n",
    "    if granules:\n",
    "        for g in granules:\n",
    "            granule_urls = ''\n",
    "            granule_poly = ''\n",
    "                       \n",
    "            # read granule title and cloud cover\n",
    "            granule_name = g['title']\n",
    "            cloud_cover = g['cloud_cover']\n",
    "    \n",
    "            # reading bounding geometries\n",
    "            if 'polygons' in g:\n",
    "                polygons= g['polygons']\n",
    "                multipolygons = []\n",
    "                for poly in polygons:\n",
    "                    i=iter(poly[0].split (\" \"))\n",
    "                    ltln = list(map(\" \".join,zip(i,i)))\n",
    "                    multipolygons.append(Polygon([[float(p.split(\" \")[1]), float(p.split(\" \")[0])] for p in ltln]))\n",
    "                granule_poly = MultiPolygon(multipolygons)\n",
    "            \n",
    "            # Get https URLs to .nc files and exclude .dmrpp files\n",
    "            granule_urls = [x['href'] for x in g['links'] if 'https' in x['href'] and '.nc' in x['href'] and '.dmrpp' not in x['href']]\n",
    "            # Add to list\n",
    "            granule_arr.append([granule_urls, cloud_cover, granule_poly])\n",
    "                           \n",
    "        page_num += 1\n",
    "    else: \n",
    "        break\n",
    " \n",
    "print(granule_arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note: At the time this tutorial was made, all 3 searches, point, bounding box, and polygon should result in the same assets being returned.\n",
    "\n",
    "### Creating a Dataframe with the resulting Links\n",
    "\n",
    "A `pandas.dataframe` can be used to store the download URLs and geometries of each file. The EMIT L2A Reflectance and Uncertainty and Mask collection contains 3 assets per granule (reflectance, reflectance uncertainty, and masks). We can see when printing this list, that there are three assets that correspond to a single polygon. For the next step we will place these into a dataframe and 'explode' the dataframe to place each of these in a separate row. If we only want a subset of these assets, we can filter them out. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# creating a pandas dataframe\n",
    "cmr_results_df = pd.DataFrame(granule_arr, columns=[\"asset_url\", \"cloud_cover\", \"granule_poly\"])\n",
    "# Drop granules with empty geometry - if any exist\n",
    "cmr_results_df = cmr_results_df[cmr_results_df['granule_poly'] != '']\n",
    "# Expand so each row contains a single url \n",
    "cmr_results_df = cmr_results_df.explode('asset_url')\n",
    "# Name each asset based on filename\n",
    "cmr_results_df.insert(0,'asset_name', cmr_results_df.asset_url.str.split('/',n=-1).str.get(-1))\n",
    "\n",
    "cmr_results_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this stage we can filter based on the assets that we want or the cloud cover. For this example lets say we are only interested in the Reflectance and the Mask. To filter by asset, we can match strings included in the asset name. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cmr_results_df = cmr_results_df[cmr_results_df.asset_name.str.contains('_RFL_') | cmr_results_df.asset_name.str.contains('MASK')]\n",
    "cmr_results_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After filtering down to the assets you want, you can output a text file with the asset urls or save the entire dataframe, then use a utility such as wget or the DAAC Data Download Tool to download the files. To download you will need to set up NASA Earthdata Login authentication using  a .netrc file. \n",
    "\n",
    "Save the asset urls to a textfile in the `/data/` folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Save text file of asset urls\n",
    "cmr_results_dfs = cmr_results_df[:-1].drop_duplicates(subset=['asset_url']) # Remove any duplicates\n",
    "cmr_results_df.to_csv('../../data/emit_asset_urls.txt', columns = ['asset_url'], index=False, header = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Downloading Files using the list of URLS/Text File\n",
    "\n",
    "To download the files using Python, you can run the cell below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define input filepath\n",
    "url_list_filepath = \"../../data/emit_asset_urls.txt\"\n",
    "# Define output directory\n",
    "output_directory = \"../../data/\"\n",
    "\n",
    "# Open Text file\n",
    "with open(url_list_filepath, \"r\") as file:\n",
    "        file_list = file.read().splitlines()\n",
    "        file.close()\n",
    "\n",
    "# EDL Authentication/Create .netrc if necessary\n",
    "earthaccess.login(persist=True)\n",
    "# Get requests https Session using Earthdata Login Info\n",
    "fs = earthaccess.get_requests_https_session()\n",
    "# Retrieve granule asset ID from URL (to maintain existing naming convention)\n",
    "for url in file_list:\n",
    "    granule_asset_id = url.split(\"/\")[-1]\n",
    "    # Define Local Filepath\n",
    "    fp = f\"{output_directory}{granule_asset_id}\"\n",
    "    # Download the Granule Asset if it doesn't exist\n",
    "    print(f\"Downloading {granule_asset_id}...\")\n",
    "    if not os.path.isfile(fp):\n",
    "        with fs.get(url, stream=True) as src:\n",
    "            with open(fp, \"wb\") as dst:\n",
    "                for chunk in src.iter_content(chunk_size=64 * 1024 * 1024):\n",
    "                    dst.write(chunk)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To download using wget, use the following in the command line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "!wget -P ../../data/ -i ../../data/emit_asset_urls.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 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: 03-22-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"
  },
  "vscode": {
   "interpreter": {
    "hash": "3292b2aceff7d39327a7519422d4180a7c9b133202090f26e797e3dd8f2c7877"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}