notebooks/measures/Search_Download_and_Plot_GeoTIFFs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1341755d",
   "metadata": {},
   "source": [
    "<center>\n",
    "    <img src='./img/nsidc_logo.png'/>\n",
    "</center>\n",
    "\n",
    "\n",
    "# **Download, Crop, Resample and Plot Multiple GeoTIFFs - Using MEaSUREs data sets as an example**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38830c9d",
   "metadata": {},
   "source": [
    "## **1. Tutorial Overview**\n",
    "\n",
    "This tutorial guides you through programmatically accessing and downloading GeoTIFF files from the NSIDC DAAC to your local computer. We then crop and resample one GeoTIFF based on the extent and pixel size of another GeoTIFF, then plot one on top of the other. \n",
    "\n",
    "We will use two data sets from the NASA [MEaSUREs](https://nsidc.org/data/measures) (Making Earth System data records for Use in Research Environments) program as an example:\n",
    "\n",
    "* [MEaSUREs Greenland Ice Mapping Project (GrIMP) Digital Elevation Model from GeoEye and WorldView Imagery, Version 2 (NSIDC-0715)](https://nsidc.org/data/nsidc-0715/versions/2)\n",
    "* [MEaSUREs Greenland Ice Velocity: Selected Glacier Site Velocity Maps from InSAR, Version 4 (NSIDC-0481)](https://nsidc.org/data/nsidc-0481/versions/4)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bb1a111",
   "metadata": {},
   "source": [
    "### **Credits**\n",
    "\n",
    "This tutorial was developed by Jennifer Roebuck of NSIDC.\n",
    "\n",
    "For questions regarding the notebook, or to report problems, please create a new issue in the [NSIDC-Data-Tutorials repo](https://github.com/nsidc/NSIDC-Data-Tutorials/issues)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f7d173fc",
   "metadata": {},
   "source": [
    "### **Objectives** \n",
    "\n",
    "1. Use the `earthaccess` library for authentication and to programmatically search for and download NSIDC DAAC data that meet specified spatial and temporal requirements. \n",
    "2. Use the `gdal` and `osr` modules from the `osgeo` package to crop and resample one GeoTIFF based on the extent and pixel size of another GeoTIFF.\n",
    "3. Use `rasterio`, `numpy` and `matplotlib` libraries to overlay one GeoTIFF on top of another.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5cee718b",
   "metadata": {},
   "source": [
    "### **Prerequisites**\n",
    "\n",
    "To run this tutorial you will need:\n",
    "1. An Earthdata Login for authentication and downloading the data. It is completely free. If you don't have one, you can register for one [here](https://urs.earthdata.nasa.gov/). \n",
    "2. A .netrc file. We recommend using a .netrc file for storing your Earthdata Login username and password. Instructions for setting one up can be found in Step 1 in this [Programmatic Data Access Guide](https://nsidc.org/data/user-resources/help-center/programmatic-data-access-guide#anchor-0). If you don't want to set one up, you will be prompted for our username and password during the tutorial.\n",
    "3. A basic understanding of python may also be helpful for this tutorial, but not necessary."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25584ec8",
   "metadata": {},
   "source": [
    "### **Example of end product**\n",
    "\n",
    "At the end of this tutorial, you will have produced a figure similar to the one below, which overlays velocity data on top of a digital elevation model:\n",
    "\n",
    "<center>\n",
    "    <img src='./img/example_geotiff_plot.png'/>\n",
    "</center>\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30fbaa69",
   "metadata": {},
   "source": [
    "### **Time requirement**\n",
    "\n",
    "This tutorial will take approximately 30 minutes to complete. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90d526d0",
   "metadata": {},
   "source": [
    "## **2. Tutorial Steps**\n",
    "\n",
    "### **Import libraries and classes**\n",
    "\n",
    "We will use the following libraries:\n",
    "1. `earthaccess` to authenticate, search and download NSIDC DAAC data \n",
    "2. `os` to list all the files we have downloaded \n",
    "3. `osgeo.gdal` and `osgeo.osr` to crop and resample one of the GeoTIFFs \n",
    "4. `rasterio` and `numpy` to read the GeoTIFFs and set up a grid for plotting the data. \n",
    "5. `matplotlib` for plotting the data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da4e53c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import earthaccess\n",
    "import os\n",
    "from osgeo import gdal, osr\n",
    "import rasterio\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from affine import Affine"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a94373e3",
   "metadata": {},
   "source": [
    "### **Authentication**\n",
    "\n",
    "We need to set up our authentication using our Earthdata Login credentials and the `login` method from `earthaccess`. If we have a .netrc with our Earthdata login credentials, we can just run the cell below and it will automatically authenticate. If we don't have a .netrc, we will be prompted to enter our Earthdata Login username and password. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "257a1e51",
   "metadata": {},
   "outputs": [],
   "source": [
    "auth = earthaccess.login()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6aad881a",
   "metadata": {},
   "source": [
    "### **Search for data using spatial and temporal filters**\n",
    "This tutorial assumes we already know which data sets we would like to download (that have data in GeoTIFF format). Each data set at NSIDC has a data set ID; it can be found after the title on the data set landing page. \n",
    "\n",
    "We will look at two data sets focused on Greenland: a Digital Elevation Model (DEM) ([NSIDC-0715](https://nsidc.org/data/nsidc-0715)) and velocity at glacier sites ([NSIDC-0481](https://nsidc.org/data/nsidc-0481)). We will use the `search_data` method from the `earthaccess` library and the following variables to search for granules within these data sets and it will output the number of granules that meet the search criteria for each data set:\n",
    "\n",
    "* `short_name` - this is the data set ID e.g., NSIDC-0715, NSIDC-0481. It can be found in the data set title on the data set landing page.\n",
    "* `version` - data set version number, also included in the data set title.\n",
    "* `cloud_hosted` - NSIDC is in the process of migrating data sets to the cloud. The data sets we are interested in are currently still archived on-prem so will set this to False.\n",
    "* `bounding_box` - sets a spatial filter by specifying latitude and longitude in the following order: W, S, E, and N.\n",
    "* `temporal` - sets a temporal filter by specifying a start and end date in the format YYYY-MM-DD.\n",
    "* `count` - this sets the maximum number of granules that will be returned in the search\n",
    "\n",
    "**NOTE** The data sets and the spatial and temporal filters used here are just examples. If there are other data sets in GeoTIFF format or different spatial and temporal filters we wish to try, then we can just modify the `short_name`, `bounding_box` and `temporal` variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62b87b2e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Search for DEM files\n",
    "print('For NSIDC-0715')\n",
    "results_dem = earthaccess.search_data(\n",
    "    short_name='NSIDC-0715',\n",
    "    version='2',\n",
    "    cloud_hosted=False,\n",
    "    bounding_box=(-33.45,68.29,-31.41,69.26),\n",
    "    temporal=('2015-12-01','2015-12-31'),\n",
    "    count=100\n",
    ")\n",
    "\n",
    "#Search for velocity data \n",
    "print('For NSIDC-0481')\n",
    "results_vel = earthaccess.search_data(\n",
    "    short_name='NSIDC-0481',\n",
    "    version='4',\n",
    "    cloud_hosted=False,\n",
    "    bounding_box=(-33.45,68.29,-31.41,69.26),\n",
    "    temporal=('2017-01-01','2017-12-31'),\n",
    "    count=100\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f58d1b42",
   "metadata": {},
   "source": [
    "### **Download the data**\n",
    "Now we have found granules that meet our search criteria we can download them using the `download` method from `earthaccess`. Note that for these particular data sets within each granule there are multiple files. So even though 1 granule was found for the DEM data set, 6 files will be downloaded as the one granule represents one tile of the DEM (tile_4_2) and for each tile there are 6 files: \n",
    "\n",
    "* _tile_4_2_30m_dem_v02.0.tif_ - the DEM\n",
    "* _tile_4_2_30m_count_v02.0.tif_ - number of individual DEMs used to estimate height at each pixel\n",
    "* _tile_4_2_30m_mindate_v02.0.tif_ - date of the most recent DEM at each pixel\n",
    "* _tile_4_2_30m_maxdate_v02.0.tif_ - date of the oldest DEM at each pixel\n",
    "* _tile_4_2_30m_mad_v02.0.tif_ - median absolute deviation (mad) of DEM heights at each pixel\n",
    "* _tile_4_2_30m_browse_v02.0.tif_ - DEM hillshade image (this is the one we will use for plotting) \n",
    "\n",
    "For the velocity data there are 20 granules, one for each time period, and there are 5 files within each granule, for example for the time period 16 July to 27 July 2017:\n",
    "\n",
    "* _TSX_E68.80N_16Jul17_27Jul17_19-41-20_vv_v04.0.tif_ - velocity magnitude \n",
    "* _TSX_E68.80N_16Jul17_27Jul17_19-41-20_vx_v04.0.tif_ - x component of velocity\n",
    "* _TSX_E68.80N_16Jul17_27Jul17_19-41-20_vy_v04.0.tif_ - y component of velocity\n",
    "* _TSX_E68.80N_16Jul17_27Jul17_19-41-20_ex_v04.0.tif_ - error of x component\n",
    "* _TSX_E68.80N_16Jul17_27Jul17_19-41-20_ey_v04.0.tif_ - error of y component \n",
    "\n",
    "Further details can be found in the data set user guides for the [DEM](https://nsidc.org/sites/default/files/documents/user-guide/nsidc-0715-v002-userguide.pdf) and [velocity](https://nsidc.org/sites/default/files/documents/user-guide/nsidc-0481-v004-userguide.pdf) data sets. \n",
    "\n",
    "The progress and size of each order will be printed out.\n",
    "\n",
    "We will create a folder called 'outputs' to store the downloaded files. This folder will be created in the folder where we are currently running this notebook. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55d7625f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Set up an outputs folder to download the data to\n",
    "path = str(os.getcwd() + '/outputs')\n",
    "if not os.path.exists(path):\n",
    "    os.mkdir(path)\n",
    "\n",
    "#Download the DEM granules \n",
    "dem_files = earthaccess.download(results_dem, path)\n",
    "\n",
    "#Download the velocity granules\n",
    "vel_files = earthaccess.download(results_vel, path)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db0008d8",
   "metadata": {},
   "source": [
    "### **Check the files that have been downloaded**\n",
    "We will list all the DEM and velocity files that were downloaded, as this is needed for the next steps. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0095777f",
   "metadata": {},
   "outputs": [],
   "source": [
    "dir_list = os.listdir(path)\n",
    "\n",
    "print('Files in outputs folder')\n",
    "\n",
    "for x in dir_list:\n",
    "    if x.endswith('.tif'):\n",
    "        print(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5afa949b",
   "metadata": {},
   "source": [
    "### **Select and read in the DEM file and velocity file**\n",
    "\n",
    "Based on the list of filenames above, we will select the GeoTIFFs that we wish to plot and input them into the cell below. We will use the 'browse' GeoTIFF of the DEM tile, as that provides the best continuous surface for visual display. For the velocity we will plot the velocity magnitude, which is denoted by 'vv' in the filename, and we will plot the velocity covering the time period 07 August to 18 August 2017. \n",
    "\n",
    "The DEM GeoTIFF covers a larger area than the velocity GeoTIFF, and we will be cropping the DEM GeoTIFF to the extent of the velocity GeoTIFF, so we will also set a filename for the cropped DEM GeoTIFF. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36027699",
   "metadata": {},
   "outputs": [],
   "source": [
    "vel_fp = str(path + '/TSX_E68.80N_07Aug17_18Aug17_19-41-22_vv_v04.0.tif')\n",
    "\n",
    "dem_fp = str(path + '/tile_4_2_30m_browse_v02.0.tif')\n",
    "\n",
    "dem_crop = str(path + '/dem_crop_100.tif')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39178989",
   "metadata": {},
   "source": [
    "### **Crop and resample DEM file based on velocity file extent and pixel size**\n",
    "\n",
    "We will use the `open` method from the `gdal` module to read the velocity GeoTIFF and get the extent and pixel size, and we will use the `SpatialReference` method from the `osr` module to get the projection information. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16184745",
   "metadata": {},
   "outputs": [],
   "source": [
    "#open the velocity GeoTIFF that we specified above\n",
    "vel_raster = gdal.Open(vel_fp)\n",
    "\n",
    "#get the grid information (extent and pixel size) from the velocity file\n",
    "geoTransform = vel_raster.GetGeoTransform()\n",
    "\n",
    "#get the projection information (EPSG code) from the velocity file\n",
    "proj=osr.SpatialReference(wkt=vel_raster.GetProjection())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80107a09-2f45-408f-90fc-85fa5923e890",
   "metadata": {},
   "source": [
    "We will now assign the information, retrieved in the cell above, to variables that can be used to crop and downsample the DEM GeoTIFF."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2db7449-d293-4d9f-a3ff-4a721c7bf6be",
   "metadata": {},
   "outputs": [],
   "source": [
    "#assign the EPSG code to a variable that will be used by GDAL to ensure the output file is in the same projection\n",
    "epsg = 'EPSG:' + proj.GetAttrValue('AUTHORITY',1)\n",
    "\n",
    "# assign the pixel size and extent info to variables that will be used to crop and downsample the data \n",
    "pixelSizeX = geoTransform[1]\n",
    "pixelSizeY = geoTransform[5]\n",
    "minx = geoTransform[0]\n",
    "maxy = geoTransform[3]\n",
    "maxx = minx + pixelSizeX * vel_raster.RasterXSize\n",
    "miny = maxy + pixelSizeY * vel_raster.RasterYSize"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "babecd8c-5ebf-4e96-ae4a-b6eff4110216",
   "metadata": {},
   "source": [
    "Now we will use these variables with the `Warp` method from the `gdal` module to crop and downsample the DEM GeoTIFF and save it as new file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2e0f586b-55d7-4eea-8193-a70c0d887ed4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#set up the keyword arguments needed for the Warp method to crop and downsample the DEM geotiff\n",
    "kwargs = { 'format': 'GTiff', 'outputBounds': [minx, miny, maxx, maxy], 'outputBoundsSRS': epsg, 'xRes': pixelSizeX, 'yRes': pixelSizeY}\n",
    "\n",
    "#use the Warp method to crop and downsample the DEM GeoTIFF and save it as a new file\n",
    "ds = gdal.Warp(dem_crop, dem_fp, **kwargs)\n",
    "ds=None"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c743d341",
   "metadata": {},
   "source": [
    "### **Set up the grids to plot the DEM and velocity data**\n",
    "To plot the cropped and downsampled DEM with the velocity data, we need to set up a grid. We will read in the data and metadata using the `open` and `meta` methods, then we will use the `transform` or `affine` methods to map the image coordinates to world coordinates, all methods from the `rasterio` package. This information will then be used to generate x and y grid locations using the `arrange` method from the `numpy` package. It will also mask the no data values from the velocity GeoTIFF to produce a nicer plot. The cell below will print out the metadata from both files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15fc1119",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read in subsetted and resampled DEM\n",
    "dem_src = rasterio.open(dem_crop)\n",
    "\n",
    "# print out metadata information\n",
    "for k in dem_src.meta:\n",
    "    print(k,dem_src.meta[k])\n",
    "\n",
    "# Retrieve the affine transformation\n",
    "if isinstance(dem_src.transform, Affine):\n",
    "     transform = dem_src.transform\n",
    "else:\n",
    "     transform = dem_src.affine\n",
    "\n",
    "N = dem_src.width\n",
    "M = dem_src.height\n",
    "dx = transform.a\n",
    "dy = transform.e\n",
    "minx = transform.c\n",
    "maxy = transform.f\n",
    "\n",
    "# Read the image data, flip upside down if necessary\n",
    "dem_crop_in = dem_src.read(1)\n",
    "if dy < 0:\n",
    "    dy = -dy\n",
    "    dem_crop_in = np.flip(dem_crop_in, 0)\n",
    "\n",
    "#Uncomment the line below if you wish to see the min/max DEM values\n",
    "#print('Data minimum, maximum = ', np.amin(data_in), np.amax(data_in))\n",
    "\n",
    "# Generate X and Y grid locations\n",
    "xdata = minx + dx/2 + dx*np.arange(N)\n",
    "ydata = maxy - dy/2 - dy*np.arange(M-1,-1,-1)\n",
    "\n",
    "dem_extent = [xdata[0], xdata[-1], ydata[0], ydata[-1]]\n",
    "\n",
    "### read in the velocity data\n",
    "vel_data = rasterio.open(vel_fp)\n",
    "\n",
    "for k in vel_data.meta:\n",
    "    print(k,vel_data.meta[k])\n",
    "\n",
    "# Retrieve the affine transformation\n",
    "if isinstance(vel_data.transform, Affine):\n",
    "     transform = vel_data.transform\n",
    "else:\n",
    "     transform = vel_data.affine\n",
    "\n",
    "N2 = vel_data.width\n",
    "M2 = vel_data.height\n",
    "dx2 = transform.a\n",
    "dy2 = transform.e\n",
    "minx2 = transform.c\n",
    "maxy2 = transform.f\n",
    "\n",
    "# Read the image data, flip upside down if necessary\n",
    "vel_in = vel_data.read(1)\n",
    "if dy2 < 0:\n",
    "    dy2 = -dy2\n",
    "    vel_in = np.flip(vel_in, 0)\n",
    "\n",
    "#Uncomment the line below if you wish to see the min/max velocity values\n",
    "#print('Data minimum, maximum = ', np.amin(vel_in), np.amax(vel_in))\n",
    "\n",
    "# Generate X and Y grid locations\n",
    "xdata2 = minx2 + dx2/2 + dx2*np.arange(N2)\n",
    "ydata2 = maxy2 - dy2/2 - dy2*np.arange(M2-1,-1,-1)\n",
    "\n",
    "vel_extent = [xdata2[0], xdata2[-1], ydata2[0], ydata2[-1]]\n",
    "\n",
    "#Need to mask the no data values in the velocity data\n",
    "vel_masked = np.ma.masked_where(vel_in == -1.0, vel_in, copy=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2febb071",
   "metadata": {},
   "source": [
    "### **Plot the DEM and velocity data**\n",
    "Now we can use the `matplotlib` library to plot the DEM with the velocity data on top. We will set the transparency of the velocity layer so we can see the DEM underneath. There is also an option to save the figure as .png; we can uncomment the last line if we want to save the image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29fa16b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "plt.figure(figsize=(8,8))\n",
    "fig = plt.imshow(dem_crop_in, extent=dem_extent, origin='lower', cmap='gray')\n",
    "fig2 = plt.imshow(vel_masked, extent=vel_extent, origin='lower', cmap='terrain', alpha=0.8)\n",
    "plt.title('Velocity and DEM')\n",
    "plt.xlabel('X (km)')\n",
    "plt.ylabel('Y (km)')\n",
    "cb = plt.colorbar(fig2, shrink=0.5)\n",
    "cb.set_label('Velocity Magnitude (m/yr)')\n",
    "\n",
    "#Option to save the figure\n",
    "#plt.savefig(\"velocity.png\", dpi=300, bbox_inches='tight', pad_inches=0.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60b01514-2442-4a51-9d50-e24eeda6550f",
   "metadata": {},
   "source": [
    "**Optional** - if you wish to clean up the outputs folder and remove all the downloaded .tif files then you can uncomment the lines in the following cell (remove the '#') and run it. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6baa3e22-12f9-4adf-b96c-ec4413ccd6d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "#for x in dir_list:\n",
    "#    if x.endswith('.tif'):\n",
    "#        os.remove(os.path.join(path, x))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ff31ea8",
   "metadata": {},
   "source": [
    "## **3. Learning Outcomes**\n",
    "\n",
    "* Search and download NSIDC DAAC data using `earthaccess`\n",
    "* Crop and resample GeoTIFF using `gdal` and `osr`\n",
    "* Overlay one GeoTIFF on another in a plot using `rasterio`, `numpy` and `matplotlib`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "182024ce",
   "metadata": {},
   "source": [
    "## **4. Additional Resources**\n",
    "\n",
    "* Further details on the `earthacess` library can be found [here](https://github.com/nsidc/earthaccess)\n",
    "* Further details on data available from NSIDC can be found [here](https://nsidc.org/data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "470d00ef",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}