BioSCape_ANG_explore.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6f23c8f0-5633-4376-b5fa-21a32e962dc7",
   "metadata": {},
   "source": [
    "## ESA 2024 Special Session - SS 08\n",
    "\n",
    "# Accessing Data from NASA’s First Biodiversity-Focused Airborne and Field Campaign, BioSCape\n",
    "**Monday, August 5, 2024, 11:45 - 1:15 PM PDT**\n",
    "- Organizer - Erin Hestir, University of California Merced\n",
    "- Co-Organizers - In Person: Anabelle Cardoso, Philip Brodrick, Michele Thornton, Erin Hestir\n",
    "- Co-Organizers - Additional: Rupesh Shrestha, Glenn Moncrieff, Adam Wilson"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74c4455b-d9e9-4955-8316-5a571a84c32e",
   "metadata": {},
   "source": [
    "## Session Overview\n",
    "\n",
    "<center><img src=\"NotebookImages/121229-87.png\"/></center>\n",
    "\n",
    "**ESA 2024 Special Session: SS 08 - Accessing Data from NASA’s First Biodiversity-Focused Airborne and Field Campaign, BioSCape**\n",
    "\n",
    "[BioSCape](https://www.bioscape.io/), the Biodiversity Survey of the Cape, is NASA’s first biodiversity-focused airborne and field campaign that was conducted in South Africa in 2023. BioSCape’s primary objective is to study the structure, function, and composition of the region’s ecosystems, and how and why they are changing. \n",
    "\n",
    "BioSCape's airborne dataset is unprecedented, with AVIRIS-NG, PRISM, and HyTES imaging spectrometers capturing spectral data across the UV, visible and infrared at high resolution and LVIS acquiring coincident full-waveform lidar. BioSCape's field dataset is equally impressive, with 18 PI-led projects collecting data ranging from the diversity and phylogeny of plants, kelp and phytoplankton, eDNA, landscape acoustics, plant traits, blue carbon accounting, and more\n",
    "\n",
    "BioSCape is committed to Open Science. All our datasets will be delivered to one of NASA's Distributed Archive and Analysis Centers. In this ESA special session we will: \n",
    "\n",
    "1) Give an overview of BioSCape’s research activities, \n",
    "2) Discuss the ways in which BioSCape data has been optimized for easier uptake by new users, and\n",
    "3) Demonstrate where and how one can access BioSCape’s field and airborne data sets.\n",
    "\n",
    "The goal of this session is to encourage wider awareness and use of BioSCape data and promote its application in research and biodiversity conservation."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a3e369c9-42df-40e3-93a2-451203c8b386",
   "metadata": {},
   "source": [
    "----\n",
    "# Tutorial: Exploring BioSCape AVIRIS-NG L2A Reflectance Data\n",
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb7828e5-d976-40f6-a805-078c94ed9bc3",
   "metadata": {},
   "source": [
    "### BioSCape AVIRIS-Next Generation Data\n",
    "\n",
    "Airborne Visible InfraRed Imaging Spectrometer - Next Generation (AVIRIS-NG) provides high signal-to-noise ratio imaging spectroscopy measurements in the solar reflected spectral range. AVIRIS-NG measures the wavelength range from 380 nm to 2510 nm with 5 nm sampling.  More information about the AVIRIS-NG instrument is [here](https://avirisng.jpl.nasa.gov/aviris-ng.html).\n",
    "\n",
    "Preliminary BioSCape AVIRIS-NG Data available here:  https://popo.jpl.nasa.gov/avng/y23_bioscape/\n",
    "\n",
    "BioSCape Airborne data interactive viewer: **BioSCape Data Portal**: https://popo.jpl.nasa.gov/mmgis-aviris/?mission=BIOSCAPE \n",
    "\n",
    "- Data are orthorectified\n",
    "- Flight lines are provided as smaller, more manageable sections of data we'll refer to as \"scenes\"\n",
    "- Within a flight line, adjacent scenes are seamless\n",
    "- To date, data are in binary/header (ENVI) formats\n",
    "- Data are calibrated to at-sensor radiance\n",
    "- `L2a Data`: Surface Reflectance and Uncertainty are available"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "906ad336-ad66-4886-b24c-c1971e94b6c0",
   "metadata": {},
   "source": [
    "### Learning Objectives\n",
    "### In this tutorial, we'll use python methods to: \n",
    "\n",
    "1. Use GDAL methods to open and explore AVIRIS-NG data\n",
    "2. Create virtual raster (vrt) to stitch adjacent scences\n",
    "3. Transform a lat/lon location to row colume locations\n",
    "4. Plot a location's spectral profile\n",
    "\n",
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97c7d120-dfd4-4d69-8d28-42c54cdfb20b",
   "metadata": {},
   "source": [
    "## AVIRIS-NG Data Files\n",
    "### **From the [BioSCape Data Portal](https://popo.jpl.nasa.gov/mmgis-aviris/?mission=BIOSCAPE), using the coordinates of the Brackenburn Reserve, we can determine and download scenes of interest**\n",
    "- Brackenburn Priviate Nature Reserve Coordinates:  **`-33.978573`**, **`23.474745`**\n",
    "\n",
    "##### Corresponding L2 Files\n",
    "- ang20231110t081307_006_L2A_OE_main_27577724_RFL_ORT_QL.tif\n",
    "- ang20231110t081307_007_L2A_OE_main_27577724_RFL_ORT_QL.tif\n",
    "- ang20231110t081307_006_L2A_OE_main_27577724_RFL_ORT.tar.gz\n",
    "- ang20231110t081307_007_L2A_OE_main_27577724_RFL_ORT.tar.gz\n",
    "\n",
    "\n",
    "| Dataset | Description | \n",
    "| -------- | --- |\n",
    "| *_RFL_ORT_QL.tif | Reflectance Quick Look Image (3 band) |\n",
    "| *-RFL_ORT | Reflectance ENVI binary file (425 band)|\n",
    "| *_RFL_ORT.hdr | Reflectance ENVI header file (txt file)|\n",
    "\n",
    "<center><img src=\"NotebookImages/ANG_naming.jpg\"/></center>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "438f5726-19c1-42ee-8acc-31eb51f167f6",
   "metadata": {},
   "source": [
    "## Study Area\n",
    "### For this tutorial, the [Brackenburn Private Nature Reserve](https://plcnetwork.co.za/member/121/Brackenburn-Private-Nature-Reserve/) was selected as a study area of interest.\n",
    "\n",
    "<center><img src=\"NotebookImages/brackenburn.PNG\"/></center>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3d5eb36-9f37-492b-9e71-72bd772dc69a",
   "metadata": {},
   "source": [
    "### Import python modules"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83871185-7114-47a0-83e6-0bf635801baf",
   "metadata": {},
   "source": [
    "- **`gdal`** a popular Geospatial Data Abstraction Library.  A translator library for raster and vector geospatial data formats.\n",
    "- **`NumPy`** is the fundamental package for scientific computing within Python and provides an N-dimensional array object suitable for multidimensional files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4a4d7b2-9148-4919-9e8f-0624d8f8879c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# python imports \n",
    "from osgeo import gdal\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from spectral.io import envi\n",
    "from pyproj import Proj\n",
    "import os\n",
    "from glob import glob"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c867bc2-0687-41bf-8a68-d37a873bbb05",
   "metadata": {},
   "source": [
    "## Examine Adjacent Quick Look geoTIFF Scenes using GDAL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af21754b-61bd-4a59-8fd1-9283c7e37246",
   "metadata": {},
   "source": [
    "### GDAL Raster Dataset\n",
    "\n",
    "- **`GDAL`** (Geospatial Data Abstraction Library) is a translator library for raster and vector geospatial data formats\n",
    "- We've already downloaded 2 adjacent scenes of AVIRIS-NG data (Quicklook and L2A Reflectance) in the Brackenburg Private Nature Reserve"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "396a0ee7-0aaf-48d5-a98d-042d5b57d96c",
   "metadata": {},
   "outputs": [],
   "source": [
    "for geotiff in glob(\"data/*.tif\"):\n",
    "    print(f\"{os.path.basename(geotiff)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ab1c80f-f1b0-4974-9409-4fc788342acc",
   "metadata": {},
   "source": [
    "### GDAL to get the number of `bands`, `rows`, and `columns` in the Quicklook files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb781d9d-5f99-4686-961d-1804e5e27293",
   "metadata": {},
   "outputs": [],
   "source": [
    "# open and read tif files for comparison\n",
    "\n",
    "ql_006 = gdal.Open('data/ang20231110t081307_006_L2A_OE_main_27577724_RFL_ORT_QL.tif')\n",
    "ql_007 = gdal.Open('data/ang20231110t081307_007_L2A_OE_main_27577724_RFL_ORT_QL.tif')\n",
    "\n",
    "nbands_006 = ql_006.RasterCount\n",
    "ncols_006 = ql_006.RasterXSize\n",
    "nrows_006 = ql_006.RasterYSize\n",
    "\n",
    "nbands_007 = ql_007.RasterCount\n",
    "ncols_007 = ql_007.RasterXSize\n",
    "nrows_007 = ql_007.RasterYSize\n",
    "\n",
    "print(f\"Bands_006:\\t{nbands_006}\")\n",
    "print(f\"Rows_006:\\t{nrows_006}\")\n",
    "print(f\"Cols_006:\\t{ncols_006}\")\n",
    "print(\"\\n\")\n",
    "print(f\"Bands_007:\\t{nbands_007}\")\n",
    "print(f\"Rows_007:\\t{nrows_007}\")\n",
    "print(f\"Cols_007:\\t{ncols_007}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "edd91a09-c275-44ba-8f85-b7ddd356ecf7",
   "metadata": {},
   "outputs": [],
   "source": [
    "ql_007.GetProjection()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f4daf1e-a60f-4644-b5e6-42bee5a9fbf8",
   "metadata": {},
   "source": [
    "### Create a Virtual Raster of the 2 Adjacent AVIRIS-NG Scenes\n",
    "\n",
    "- A **virtual raster** (VRT) is a GDAL method of combining multiple raster tiles into one file. It enables visualisation of whole raster datasets and faster navigation around your dataset in GIS software.\n",
    "- Using 2 adjacent flight line AVIRIS-NG quicklook images (3-band true color geoTIFF file) we'll use GDAL to create a virtual raster (vrt). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc421149-1a49-43a7-bb07-e62bb6e89299",
   "metadata": {},
   "outputs": [],
   "source": [
    "#build vrt - using AVIRIS-NG quicklook images \n",
    "\n",
    "# running gdalbuilvrt as a gdal command line using the '!' in front says do this as a command, not as python\n",
    "# using scences from the Bracnenburn Private Nature Reserve\n",
    "\n",
    "!gdalbuildvrt brack_mosaic.vrt data/*.tif\n",
    "\n",
    "# read in vrt, ql ~ quicklook\n",
    "brack_vrt = gdal.Open('brack_mosaic.vrt')\n",
    "\n",
    "nbands_vrt = brack_vrt.RasterCount\n",
    "ncols_vrt = brack_vrt.RasterXSize\n",
    "nrows_vrt = brack_vrt.RasterYSize\n",
    "\n",
    "print(f\"Bands_vrt:\\t{nbands_vrt}\")\n",
    "print(f\"Rows_vrt:\\t{nrows_vrt}\")\n",
    "print(f\"Cols_vrt:\\t{ncols_vrt}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64cca1db-9be1-439b-bdd2-ae231f7166ea",
   "metadata": {},
   "source": [
    "## Examine AVIRIS-NG Quiklook Reflectance Data as a Numerical Array\n",
    "\n",
    "GDAL is a powerful library when it comes to accessing geospatial raster data, but it does not provide many functionalities for doing calculations. For more advanced computing, we will **read in the raster data** as a numerical array in order to use the capabilities in the NumPy library.\n",
    "\n",
    "You can convert an existing Gdal Dataset or a Band into a **numpy array** with the **`ReadAsArray()`** function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d40d627e-5bfa-44ed-ae07-8eb98437f5de",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Recall that we created the quicklook virtual raster (ql_vrt) from two adjacent scenes\n",
    "brack_array = brack_vrt.ReadAsArray()\n",
    "\n",
    "# look at the shape. GDAL reads the array as band,x,y.  \n",
    "print(brack_array.shape)                     "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "674e9b30-0f47-4334-83ef-0250b74fc5ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "#To plot the array with matplotlib, the data need to be arranged as y, x, band.\n",
    "# We'll transpose the data\n",
    "trans_array = brack_array.transpose((1,2,0))\n",
    "print(trans_array.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49dd3216-2567-4e27-a0de-bfd9da6ef273",
   "metadata": {},
   "outputs": [],
   "source": [
    "# To demonstrate that data needs transposed, this will give an error because the data needs to be transposed\n",
    "# This is the error:  TypeError: Invalid shape (3, 1268, 1216) for image data\n",
    "#plt.imshow(brack_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4542de6-1d75-491b-999f-33c34076f8c2",
   "metadata": {},
   "source": [
    "### The imshow() function in the pyplot module of matplotlib library is used to display data as an image. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b96d4718-1d2b-442a-b3f4-d30f679fd01f",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = [10,7]\n",
    "plt.imshow(trans_array)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "393aaf8c-52be-4e1c-a0a6-5270bdd4b42b",
   "metadata": {},
   "source": [
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de36f05b-4f63-41df-87ec-5ddd14520476",
   "metadata": {},
   "source": [
    "##  AVIRIS-NG L2A Reflectance DataCube - Image Spectroscopy (ENVI format)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e074c2c2-c6ec-4304-81f3-109deedde1c1",
   "metadata": {},
   "source": [
    "#### Recall that we examined 2 adjacent tiles from the same flight path for the following files\n",
    "- Recall we downloaded files \n",
    "- ang20231110t081307_006_L2A_OE_main_27577724_RFL_ORT_QL.tif\n",
    "- ang20231110t081307_007_L2A_OE_main_27577724_RFL_ORT_QL.tif\n",
    "\n",
    "##### Corresponding Reflectance Files\n",
    "- ang20231110t081307_006_L2A_OE_main_27577724_RFL_ORT.tar.gz\n",
    "- ang20231110t081307_007_L2A_OE_main_27577724_RFL_ORT.tar.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65f7cbfa-a618-46f8-8058-30bfb78c3eac",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Need to untar/uncompress data files\n",
    "#!tar -zxvf aang20231110t081307_006_L2A_OE_main_27577724_RFL_ORT.tar.gz\n",
    "#!tar -zxvf ang20231110t081307_007_L2A_OE_main_27577724_RFL_ORT.tar.gz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8ad1ba6e-9427-4ac1-9c85-77b9426631b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "for files in glob(\"data/*\"):\n",
    "    print(f\"{os.path.basename(files)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7eea72cd-c325-483d-acea-b0e30a087a14",
   "metadata": {},
   "source": [
    "----\n",
    "### The current AVIRIS-NG Reflectance Files are in ENVI file formats which are binary/header file pairs\n",
    "\n",
    "### Let's examine the header file of one of the ENVI scenes that we downloaded:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fab5f00-1701-4fdb-aeca-0e448fb9cb15",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "hdr_f = 'data/ang20231110t081307_007_L2A_OE_main_27577724_RFL_ORT.hdr'\n",
    "with open(hdr_f, mode='r') as hdr:\n",
    "    lines = (hdr.read())\n",
    "    print(lines)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f678d644-ed24-43cd-ab53-e68525a9aae4",
   "metadata": {},
   "source": [
    "### Open the ENVI Reflectance File as a GDAL raster dataset and print the dimensions:\n",
    "#### GDAL to get the number of bands, rows, and columns in the Reflectance file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b49e93d8-165a-4def-9738-04e7ede5a370",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open the ENVI file and read the file bands, row, cols\n",
    "\n",
    "rfl_007_open = gdal.Open('data/ang20231110t081307_007_L2A_OE_main_27577724_RFL_ORT')\n",
    "\n",
    "nbands = rfl_007_open.RasterCount\n",
    "nrows = rfl_007_open.RasterYSize\n",
    "ncols = rfl_007_open.RasterXSize\n",
    "\n",
    "print(f\"Bands:\\t{nbands}\")\n",
    "print(f\"Rows:\\t{nrows}\")\n",
    "print(f\"Cols:\\t{ncols}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1df1d90e-dd94-48cb-80c7-644182aace2a",
   "metadata": {},
   "source": [
    "#### GDAL to get Metadata and Projection Information\n",
    "- The GDAL **GetMetadata** function can read some metadata from images and organizes the metadata in the form of a dictionary\n",
    "- The GDAL **GetProjection** function to obtain file projection information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a180a7fb-7fd5-40f2-9700-6d553f7d33c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "rfl_007_open.GetMetadata()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "17f9f1cb-2937-40bb-a972-8a8ad6fa2381",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"ENVI image WKT: \\n\"+str(rfl_007_open.GetProjection()))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c85df563-b834-46fc-9237-bab8c8f7b773",
   "metadata": {},
   "source": [
    "### UTM zone 34S, datum WGS-84 = EPSG Code: 32734 "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c363fcc-442e-4cad-a9c8-b708417b222e",
   "metadata": {},
   "source": [
    "#### Read one spectral band array and plot it\n",
    "- Red:  Band 57 is ~657nm (the center of the Landsat red band)\n",
    "- Green: Band 38 is ~562nm (the center of the Landsat green band)\n",
    "- Blue:  Band 22 is ~482nm (the center of the Landsat blue band)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bbaf2b9e-3ff1-4dc5-b172-c65b3be9b916",
   "metadata": {},
   "outputs": [],
   "source": [
    "#### uncomment one of the GetRaseterBand lines to run\n",
    "img_red = rfl_007_open.GetRasterBand(57).ReadAsArray()  # Band 57 is 657nm (the center of the Landsat red band)\n",
    "#img_green = rfl_007_open.GetRasterBand(38).ReadAsArray()  # Band 38 is 562nm (the center of the Landsat green band)\n",
    "#img_blue = rfl_007_open.GetRasterBand(22).ReadAsArray()  # Band 22 is 482nm (the center of the Landsat blue band)\n",
    "plt.rcParams['figure.figsize'] = [7,17]\n",
    "plt.rcParams['figure.dpi'] = 100\n",
    "plt.imshow(img_red, vmin=0, vmax=.20)\n",
    "#plt.imshow(img_green, vmin=0, vmax=0.20)\n",
    "#plt.imshow(img_blue, vmin=0, vmax=0.20)\n",
    "plt.colorbar(fraction=0.04, pad=0.04)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31f97124-f575-4ed2-b414-80573b25150d",
   "metadata": {},
   "source": [
    "## Create a Red, Green, Blue Composite and Visualize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "906a91ea-2d57-41fa-a002-374b870c0f72",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# get r,g,b arrays\n",
    "red, green, blue = rfl_007_open.GetRasterBand(57).ReadAsArray(), rfl_007_open.GetRasterBand(38).ReadAsArray(), rfl_007_open.GetRasterBand(22).ReadAsArray()\n",
    "\n",
    "# set fill values (-9999.) to 0 for each array\n",
    "red[red == -9999.], green[green == -9999.], blue[blue == -9999.] = 0, 0, 0\n",
    "\n",
    "# function scales reflectance values to 8 bits\n",
    "scale8bit = lambda a: ((a - a.min()) * (1/(a.max() - a.min()) * 255)).astype('uint8')\n",
    "\n",
    "# get 8bit arrays for each band\n",
    "red8, green8, blue8 = scale8bit(red), scale8bit(green), scale8bit(blue)\n",
    "\n",
    "# set rescaled fill pixels back to 0 for each array\n",
    "red8[red == 0], green8[green == 0], blue8[blue == 0] = 0, 0, 0\n",
    "\n",
    "rgb_stack = np.zeros((nrows,ncols,3),'uint8')\n",
    "rgb_stack[...,0], rgb_stack[...,1], rgb_stack[...,2] = red8, green8, blue8\n",
    "\n",
    "plt.rcParams['figure.figsize'] = [7,17]\n",
    "plt.rcParams['figure.dpi'] = 100\n",
    "plt.imshow(rgb_stack, vmin=0, vmax=.40)\n",
    "plt.colorbar(fraction=0.04, pad=0.04)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2230b543-600a-44b4-90ef-9e29574d1f57",
   "metadata": {},
   "source": [
    "## Apply a Histogram Stretch to for an Improved Visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "51324bd3-e60c-4ea5-96a8-4bb73d7ac886",
   "metadata": {},
   "outputs": [],
   "source": [
    "rgb_stack.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "307827d6-ca24-4933-9346-566bc1b72879",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# apply histogram equalization to each band\n",
    "for i in range(rgb_stack.shape[2]):\n",
    "\n",
    "    # band i\n",
    "    b = rgb_stack[:,:,i]\n",
    "    \n",
    "    # histogram from flattened (1d) image\n",
    "    b_histogram, bins = np.histogram(b.flatten(), 256)\n",
    "\n",
    "    # cumulative distribution function\n",
    "    b_cumdistfunc = b_histogram.cumsum()\n",
    "\n",
    "    # normalize\n",
    "    b_cumdistfunc = 255 * b_cumdistfunc / b_cumdistfunc[-1]\n",
    "\n",
    "    # get new values by linear interpolation of cdf\n",
    "    b_equalized = np.interp(b.flatten(), bins[:-1], b_cumdistfunc)\n",
    "    \n",
    "    # reshape to 2d and add back to rgb_stack\n",
    "    rgb_stack[:,:,i] = b_equalized.reshape(b.shape)\n",
    "\n",
    "plt.rcParams['figure.figsize'] = [7,17]\n",
    "plt.rcParams['figure.dpi'] = 100\n",
    "plt.imshow(rgb_stack, vmin=0, vmax=.40)\n",
    "plt.colorbar(fraction=0.04, pad=0.04)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39120b05-d473-45dd-a729-672f7882f8ca",
   "metadata": {},
   "source": [
    "-----\n",
    "-----\n",
    "\n",
    "## Plot a Spectral Profile From a Point of Interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a852f15-dece-4c13-9294-eebe586b93d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from spectral.io import envi\n",
    "import numpy as np\n",
    "from pyproj import Proj\n",
    "import os\n",
    "from osgeo import gdal"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13f713b7-ab1e-4cb6-9366-a8eb96bd1431",
   "metadata": {},
   "source": [
    "### We'll plot the spectral profiles for a pixel from the Brackenburg Reserve.\n",
    "\n",
    "Recall the Brackenburg Private Nature Reserve geographic location provided on their website\n",
    "\n",
    "**[-33.978573 lat, 23.474745 lon]**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b48e3dc5-1750-4011-b7e0-30a7c74f7443",
   "metadata": {},
   "outputs": [],
   "source": [
    "# let's define the coordinates and data file from which we'll extract spectra and the lat/lon coordinates\n",
    "data_file = 'data/ang20231110t081307_007_L2A_OE_main_27577724_RFL_ORT'\n",
    "coords = [-33.978573, 23.474745]\n",
    "coords"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19ca2077-38ab-4733-8903-576fd0bd2445",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_file"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff1980a9-fdba-4dbf-9e4f-5aee63891578",
   "metadata": {},
   "source": [
    "## Geotransform\n",
    "\n",
    "## Before we can extract pixel spectral information, we need to convert the lat/lon coordinate to an image coordinate (pixel, line) space\n",
    "\n",
    "### - The point location is in a geographic coordinate system\n",
    "### - We saw earlier, the image file is in a UTM zone 34S coordinate reference system\n",
    "\n",
    "- 1. project the lat/lon point from geographic (units = decimal degrees) to the images projected (x,y) UTM coordinates (units = meters)\n",
    "- 2. projected x,y (meter) coordinates need to be translated into a gridded offset of pixel, line values\n",
    "- 3. extract pixel band values for the x,y offset a position \n",
    "\n",
    "#### 1. Project point location lat/lon to UTM (zone 34S): Recall, the AVIRIS-NG Scene is in UTM zone 34S (EPSG:32734)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0162ceca-73ab-4618-bf53-b604703331ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "# step to translate coordinates\n",
    "# EPSG:32734 - WGS-84 / UTM zone 34S\n",
    "\n",
    "from pyproj import Proj\n",
    "p = Proj(\"EPSG:32734\", preserve_units=False)\n",
    "x,y = p(23.474745, -33.978573)\n",
    "print('x value in UTM meters:', x)\n",
    "print('y value in UTM meters:', y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f92e9a8e-f51c-446e-a090-993206590062",
   "metadata": {},
   "outputs": [],
   "source": [
    "coords_UTM = [6237459.364997024, 728621.594343807]\n",
    "\n",
    "# define aonther location for comparison\n",
    "coords_UTM2 = [6239776.9, 730509.6]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "28d5c241-9250-43c9-b1cf-24e5b0418a5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "gdal_ds = gdal.Open(data_file)\n",
    "proj_native = gdal_ds.GetProjection()\n",
    "proj_native"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ac4cabf2-d829-43b8-82b0-29f14cbff19d",
   "metadata": {},
   "source": [
    "#### 2. And now let's get and print AVIRIS-NG Scene information that we need in order to transform the UTM point x,y values to an image row/col location"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "919514fa-61c0-4f64-a774-16d53b90354a",
   "metadata": {},
   "outputs": [],
   "source": [
    "trans = gdal_ds.GetGeoTransform() # x_ul, x_px, x_rot, y_ul, y_rot, y_px\n",
    "print('\\nGetGeoTransform Return Values =', trans) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1349da94-137f-4424-98bd-4214eb543b41",
   "metadata": {},
   "source": [
    "#### GetGeotransform\n",
    "**GT** is the geotransform acquired with img.GetGeoTransform()\n",
    "\n",
    "-  GT(0): X origin, x_ul, **`730461.8657409984`**\n",
    "-  GT(1): X resolution in the pixel space, x_px, **`5.0`**\n",
    "-  GT(2). Represent the rotation of the pixel space from the geodetic space, x_rot, **`0.0`**\n",
    "-  GT(3). Y origin, y_ul, **`6240060.191764411`**\n",
    "-  GT(4). Represent the rotation of the pixel space from the geodetic space, y_rot, **`0.0`**\n",
    "-  GT(5): -1 * Y resolution in the pixel space, y_px, **`-5.0`**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8928f4f-0cea-49d1-b2d4-6d72194fd2f7",
   "metadata": {},
   "source": [
    "#### 3. x,y offset values are a position on the x,y grid of the pixel of interest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b28c602d-a918-4f6c-a14a-4250f86dcc21",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_px_offset = int(round((coords_UTM[1] - trans[0]) / trans[1]))\n",
    "y_px_offset = int(round((coords_UTM[0] - trans[3]) / trans[5]))\n",
    "\n",
    "print(x_px_offset, y_px_offset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b0eedc9d-7d19-4c96-aa58-ce2a91c1bb0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "x_px_offset2 = int(round((coords_UTM2[1] - trans[0]) / trans[1]))\n",
    "y_px_offset2 = int(round((coords_UTM2[0] - trans[3]) / trans[5]))\n",
    "\n",
    "print(x_px_offset2, y_px_offset2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ea7b292-4e21-4eaf-b3cf-eacbd619be0e",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = envi.open(os.path.splitext(data_file)[0] + '.hdr')\n",
    "ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db8c5994-817a-4fc6-9c16-57cc029b0e9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "offset_size = 1\n",
    "data = ds.open_memmap(interleave='bip')[y_px_offset - 1: y_px_offset + 1, x_px_offset - 1: x_px_offset + 1, :]\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12b9017e-f0e6-4b04-9459-da88af89b284",
   "metadata": {},
   "outputs": [],
   "source": [
    "offset_size = 1\n",
    "data2 = ds.open_memmap(interleave='bip')[y_px_offset2 - 1: y_px_offset2 + 1, x_px_offset2 - 1: x_px_offset2 + 1, :]\n",
    "data2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94d3f2fc-0f84-4459-9c1e-4cb3d174feee",
   "metadata": {},
   "outputs": [],
   "source": [
    "wl = np.array([float(x) for x in ds.metadata['wavelength']])\n",
    "wl2 = np.array([float(x) for x in ds.metadata['wavelength']])\n",
    "\n",
    "plt.plot(wl, np.mean(data, axis=(0, 1)))\n",
    "plt.rcParams['figure.figsize'] = [8, 8]\n",
    "plt.xlabel('Wavelength', fontsize=20)\n",
    "plt.ylabel('Reflectance', fontsize=20)\n",
    "plt.show()\n",
    "\n",
    "# plot the average spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecce438d-4748-4a5a-b919-f575da8cc6d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a list of wavelengths that are \"bad\" \n",
    "\n",
    "bblist = np.ones((425,))  # create a 1D array with values ones\n",
    "# set tails and atmospheric window to zero\n",
    "bblist[0:14] = 0        # tail\n",
    "bblist[189:225] = 0     # atmospheric window\n",
    "bblist[281:336] = 0     # atmospheric window\n",
    "bblist[405:] = 0        # tail"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "349c793e-a8ee-4ceb-9ac2-5da860860227",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.rcParams['figure.figsize'] = [8, 8]\n",
    "\n",
    "wl[bblist == 0] = np.nan \n",
    "wl2[bblist == 0] = np.nan\n",
    "plt.plot(wl, np.mean(data, axis=(0, 1)), color = 'g')\n",
    "plt.plot(wl2, np.mean(data2, axis=(0, 1)), color = 'r')\n",
    "plt.xlabel('Wavelength', fontsize=20)\n",
    "plt.ylabel('Reflectance', fontsize=20)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97fb5ddc-6346-4ad2-b1ca-5290232ecabc",
   "metadata": {},
   "source": [
    "## Congrats! - You're on your way to exploring AVIRIS-NG data from the BioSCape Project"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f494b587-4b88-4813-8e3d-c272851321ef",
   "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}