python/how-tos/How_to_Orthorectify.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How To: Orthorectify EMIT Data\n",
    "\n",
    ">**Warning: This notebook uses a lot of memory, if using 2i2c, you will need a 29.7 GB, 3.75 CPU instance.**\n",
    "\n",
    "EMIT Data is provided in non-orthorectified format to reduce data size. The `location` group contains latitude and longitude values of each pixel as well as a geometric lookup table (GLT) that can be used to orthocorrect the imagery. The GLT is an array that provides relative downtrack and crosstrack reference locations from the raw scene to facilitate fast projection of the dataset.\n",
    "\n",
    "This notebook will walk through using the `emit_tools` module to do the orthorectification process as well as a in depth walkthrough of the process.\n",
    "\n",
    "If you are planning to convert the EMIT `netCDF4` files to `.envi` format, the `reformat.py` available in the [emit-sds/emit-utils](https://github.com/emit-sds/emit-utils) repository can be used to orthorectify the imagery during the conversion process.\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 open an EMIT file as an `xarray.Dataset`\n",
    "+ How to orthorectify an EMIT file"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Using emit_xarray to Orthorectify\n",
    "\n",
    "Import the required Python libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Packages\n",
    "import os\n",
    "import earthaccess\n",
    "import netCDF4 as nc\n",
    "from osgeo import gdal\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import hvplot.xarray\n",
    "import holoviews as hv\n",
    "import sys\n",
    "sys.path.append('../modules/')\n",
    "from emit_tools import emit_xarray"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Login to your NASA Earthdata account and create a `.netrc` file using the `login` function from the `earthaccess` library. If you do not have an Earthdata Account, you can create one [here](https://urs.earthdata.nasa.gov/home). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "earthaccess.login(persist=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this notebook we will download the files necessary using `earthaccess`. You can also access the data in place or stream it, but this can slow due to the file sizes. Provide a URL for an EMIT L2A Reflectance granule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "url = 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/EMITL2ARFL.001/EMIT_L2A_RFL_001_20220903T163129_2224611_012/EMIT_L2A_RFL_001_20220903T163129_2224611_012.nc'"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get an HTTPS Session using your earthdata login, set a local path to save the file, and download the granule asset - This may take a while, the reflectance file is approximately 1.8 GB. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 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",
    "granule_asset_id = url.split('/')[-1]\n",
    "# Define Local Filepath\n",
    "fp = f'../../data/{granule_asset_id}'\n",
    "# Download the Granule Asset if it doesn't exist\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)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the `emit_xarray` function imported from the `emit_tools` module to read in the data.\n",
    "\n",
    "The `emit_xarray` function will open and and place data from the groups of an EMIT dataset into an `xarray.Dataset`. It includes an `ortho` option which is set to `False` by default. This function reads the root group, `location` group, and `sensor_band_parameters` groups from the EMIT dataset. It then flattens all of the variables contained within those 3 groups into a single `xarray_dataset`. When `ortho=True`, it uses the GLT layers to build a lat/lon grid and places the values from the `reflectance` root group into the grid. The lat and lon grid is used as dimensional coordinates in the output `xarray.dataset` while the `wavelengths` and `fwhm` from the `sensor_band_parameters` group are used as non-dimensional coordinates. The `reflectance` data from the root group of the netCDF is then orthorectified and added as a variable in the `xarray.Dataset`.\n",
    "\n",
    "Read in a dataset using the `emit_xarray` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = emit_xarray(fp, ortho=True)\n",
    "ds"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize a red band. First assign fill-values to `np.nan`, then select and plot a band using the `sel` function from `xarray` and the `hvplot` library. When we orthorectify (place data on a grid and rotate it) we end up with a bunch of extra pixels that have no data. These are filled with -9999, which will distort our visualization unless we assign them to `np.nan` or mask them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assign fill-values to nan\n",
    "ds.reflectance.data[ds.reflectance.data == -9999] = np.nan\n",
    "# Get red band and plot example\n",
    "ds.sel(wavelengths=650,method='nearest').hvplot.image(cmap='viridis', aspect = 'equal', frame_width=500, rasterize=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The orthorectified dataset can also be saved as a netCDF4 output that can be reopened using the `xarray.open_dataset` function if so desired using the cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.to_netcdf('../../data/example_emit_ortho_out.nc')\n",
    "# Example for Opening \n",
    "# ds = xr.open_dataset('../data/example_emit_ortho_out.nc')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now open a dataset with `ortho=False`. Notice that the data is not on a Lat/Lon grid, instead each downtrack and crosstrack pixel center is geotagged with a latitude and longitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = emit_xarray(fp,ortho=False)\n",
    "ds"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see the difference in structure and the plotted data if the dataset is orthorectified or not. It's also not necessary to mask fill-values generated by the orthorectification process."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get red band and plot example\n",
    "ds.sel(wavelengths=850,method='nearest').hvplot.image(cmap='viridis', aspect = 'equal', frame_width=500, rasterize=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Detailed Explanation of Orthorectification and xarray.Dataset Construction\n",
    "\n",
    "As mentioned, the orthorectification process for EMIT data utilizes a premade geometric lookup table (GLT) that is included in the `location` group of the file. `glt_x` and `glt_y` contain the pixel indices from the raw dataset to construct an orthocorrected array. These indices data can be pulled from the raw dataset to populate an orthocorrected array using numpy broadcasting.\n",
    "\n",
    "To do this, first read in the `location` group data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loc = xr.open_dataset(fp, engine = 'h5netcdf', group='location')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, build a numpy array using the GLT arrays from the `location` group. Use the `stack` function to stack the x and y GLT arrays and use the function `nan_to_num` to set the NaN values to the `GLT_NODATA_VALUE`, which is 0 for EMIT datasets. After doing this, you can check the shape of the array to see the expected height and width in pixels of the image after orthorectification."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "GLT_NODATA_VALUE=0\n",
    "glt_array = np.nan_to_num(np.stack([loc['glt_x'].data,loc['glt_y'].data],axis=-1),nan=GLT_NODATA_VALUE).astype(int)\n",
    "glt_array.shape"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After we have a GLT array, we need to open and create a numpy array from the reflectance dataset. Open the dataset using `xarray` and assign the data to an array. Here you can check the shape of the array to see the dimensions of the uncorrected data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = xr.open_dataset(fp,engine = 'h5netcdf')\n",
    "ds_array = ds['reflectance'].data\n",
    "ds_array.shape"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have both arrays, we need to assign a fill value for positions that have no data. We can then construct an empty array with the dimensions we desire (GLT shape) and populate it with values from the dataset. To do this we first use np.zeros to create an array of all zeros then add the `fill_value` to the array of zeros to change the zeros to the `fill_value`.\n",
    "\n",
    "Next we can build an array of valid locations by omitting the portions of the array containing the `GLT_NODATA_VALUES`.\n",
    "\n",
    "After that, change the indexing of the GLT array which is one based to zero based by subtracting 1 from them.\n",
    "\n",
    "Lastly, we can use broadcasting/indexing to pull through the values we want from the dataset into our new output array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Build Output Dataset\n",
    "fill_value = -9999\n",
    "out_ds = np.zeros((glt_array.shape[0], glt_array.shape[1], ds_array.shape[-1]), dtype=np.float32) + fill_value\n",
    "valid_glt = np.all(glt_array != GLT_NODATA_VALUE, axis=-1)\n",
    "# Adjust for One based Index\n",
    "glt_array[valid_glt] -= 1 \n",
    "# Use indexing/broadcasting to populate array cells with 0 values\n",
    "out_ds[valid_glt, :] = ds_array[glt_array[valid_glt, 1], glt_array[valid_glt, 0], :]\n",
    "out_ds.shape"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After this step, we have an array of values that is orthorectified, but if we want to have a grid of lat/lon values we still need to calculate them using the geotransform. We can retrieve this information from the dataset attributes (`ds.attrs['geotransform']`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "GT = ds.geotransform\n",
    "GT"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The geotransform consists of 2 formulas with 6 coefficients used to calculate the xy-grid to project the dataset. \n",
    "\n",
    "`x_geo = GT[0] + x * GT[1] + y * GT[2]`  \n",
    "\n",
    "`y_geo = GT[3] + x * GT[4] + y * GT[5]`\n",
    "\n",
    "`GT[0]`  - The x-coordinate of the upper-left corner of the upper-left pixel  \n",
    "`GT[1]`  - W-E pixel width  \n",
    "`GT[2]`  - Row rotation (zero for North up images)  \n",
    "`GT[3]`  - The y-coordinate of the upper-left corner of the upper-left pixel  \n",
    "`GT[4]`  - Column rotation (zero for North up images)  \n",
    "`GT[5]`  - N-S pixel height (negative value for a north-up image)  \n",
    "\n",
    "Create empty arrays for the x and y (lon and lat) based on the dimensions of the dataset.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create Array for Lat and Lon and fill\n",
    "dim_x = loc.glt_x.shape[1]\n",
    "dim_y = loc.glt_x.shape[0]\n",
    "lon = np.zeros(dim_x)\n",
    "lat = np.zeros(dim_y)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After we have an empty array with the correct dimensions, we can populate it with values using the geotransform formula to build a lat/lon grid for plotting the orthorectified data. We can remove the `GT[2]` and `GT[4]` terms from the equation since our orthorectified image is North up. We also want to shift the geotransform by half a pixel since we want to build an array of pixel centers. We do this by adding half a pixel width to `GT[0]` and `GT[3]`.  Calclutate the latitude and longitude for each row (x_dim) and column (y_dim) of the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for x in np.arange(dim_x):\n",
    "    x_geo = (GT[0]+0.5*GT[1]) + x * GT[1]\n",
    "    lon[x] = x_geo\n",
    "for y in np.arange(dim_y):\n",
    "    y_geo = (GT[3]+0.5*GT[5]) + y * GT[5]\n",
    "    lat[y] = y_geo\n",
    "\n",
    "lat, lon"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a last step we can place this and remaining information from the EMIT dataset, such as attributes and wavelengths into an `xarray.Dataset` if desired. To do this, build dictionaries for the coordinates and variables components of the dataset. In this example, the variable would be `reflectance`. First, we'll want to read in the `sensor_band_parameters` group to get the wavelength and fwhm info, then build the dictionaries.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wvl = xr.open_dataset(fp, engine = 'h5netcdf', group='sensor_band_parameters')\n",
    "coords = {'lat':(['lat'],lat), 'lon':(['lon'],lon), **wvl.variables} ## ** upacks the existing dictionary from the wvl dataset.\n",
    "data_vars = {'reflectance':(['lat','lon','bands'], out_ds)}"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now build an `xarray.Dataset` using the dictionaries, then add attributes to the dataset from the non-orthocorrected data. We also can add the crs using the `rio.write_crs` function, to allow `rasterio` and `rioxarray` to recognize the CRS for future use.\n",
    "\n",
    ">Note: the wavelength and fwhm are included as non-dimensional coordinates because they are independent of Lat/Lon."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_xr = xr.Dataset(data_vars=data_vars, coords=coords, attrs= ds.attrs)\n",
    "out_xr['reflectance'].attrs = ds['reflectance'].attrs\n",
    "out_xr.coords['lat'].attrs = loc['lat'].attrs\n",
    "out_xr.coords['lon'].attrs = loc['lon'].attrs\n",
    "out_xr.rio.write_crs(ds.spatial_ref,inplace=True) # Add CRS in easily recognizable format\n",
    "out_xr"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have an orthorectified dataset with coordinates that can be plotted. We can also mask the `fill_values` of -9999 to make a clearer figure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_xr['reflectance'].data[out_xr['reflectance'].data == fill_value] = np.nan\n",
    "out_xr.sel(bands=63).hvplot.image(cmap='viridis', aspect = 'equal', frame_width=500, rasterize=True)"
   ]
  },
  {
   "attachments": {},
   "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: 05-30-2023  \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
}