python/tutorials/Working_with_EMIT_L2B_Mineralogy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "28e9239d-5a33-4e90-b1e8-d18bc44f4f45",
   "metadata": {},
   "source": [
    "> **Disclaimer:** The [**EMITL2BMIN**](https://doi.org/10.5067/EMIT/EMITL2BMIN.001) product is generated to support the EMIT mission objectives of constraining the sign of dust related radiative forcing. Ten mineral types are the core focus of this work: Calcite, Chlorite, Dolomite, Goethite, Gypsum, Hematite, Illite+Muscovite, Kaolinite, Montmorillonite, and Vermiculite.  The [**EMIT_L3_ASA**]() product contain the aggregate abundance of these minerals at a coarser resolution for use in Earth System Models. Additional minerals are included in the **EMITL2BMIN** product for transparency but were not the focus of this product. Further validation is required to use these additional mineral maps, particularly in the case of resource exploration. Similarly, the separation of minerals with similar spectral features, such as a fine-grained goethite and hematite, is an area of active research. The results presented here are an initial offering, but the precise categorization is likely to evolve over time, and the limits of what can and cannot be separated on the global scale is still being explored. The user is encouraged to read the [Algorithm Theoretical Basis Document (ATBD)](https://lpdaac.usgs.gov/documents/1659/EMITL2B_ATBD_v1.pdf) for more details.\n",
    "\n",
    "# Working with EMIT L2B Mineralogy Data\n",
    "\n",
    "**Summary**  \n",
    "\n",
    "In this notebook we will open the [EMIT L2B Estimated Mineral Identification and Band Depth and Uncertainty (EMITL2BMIN)](https://doi.org/10.5067/EMIT/EMITL2BMIN.001) products, find a mineral of interest from the ten mineral types focused on by the EMIT mission, evaluate uncertainty, orthorectify the data, then create an output mask or vector file for the granule.\n",
    "\n",
    "**Background**\n",
    "\n",
    "The EMIT instrument is an imaging spectrometer that measures light in visible and infrared wavelengths. These measurements display unique spectral signatures that correspond to the composition on the Earth's surface. The EMIT mission focuses specifically on mapping the composition of minerals to better understand the effects of mineral dust throughout the Earth system and human populations now and in the future. More details about EMIT and its associated products can be found in the **README.md** and on the [EMIT website](https://earth.jpl.nasa.gov/emit/).\n",
    "\n",
    "The [EMITL2BMIN](https://doi.org/10.5067/EMIT/EMITL2BMIN.001) data product provides estimated mineral identification, band depths and uncertainty in a spatially raw, non-orthocorrected format. Two spectral groups, which correspond to different regions of the spectra, are identified independently and often co-occur are used to identify minerals. These estimates are generated using the [Tetracorder system](https://www.usgs.gov/publications/tetracorder-user-guide-version-44?_gl=1*1eoj33d*_ga*MTU3MTA3ODgxNS4xNjQ5MTg1MDgx*_ga_0YWDZEJ295*MTY4NjkyNTg0Mi40NC4xLjE2ODY5MjU4NzMuMC4wLjA.)([code](https://github.com/PSI-edu/spectroscopy-tetracorder)) and are based on [EMITL2ARFL](https://doi.org/10.5067/EMIT/EMITL2ARFL.001) reflectance values. The product also consists of an EMIT_L2B_MINUNCERT file, which provides band depth uncertainty estimates calculated using surface Reflectance Uncertainty values from the [EMITL2ARFL](https://doi.org/10.5067/EMIT/EMITL2ARFL.001) data product. The band depth uncertainties are presented as standard deviations, and the fit score for each mineral identification is also provided as the coefficient of determination (r<sup>2</sup>) of the match between the continuum normalized library reference and the continuum normalized observed spectrum. Associated metadata indicates the name and reference information for each identified mineral, and additional information about aggregating minerals into different categories, and the code used for product generation is available in the [emit-sds-l2b repository]().\n",
    "\n",
    "**Disclaimer**\n",
    "\n",
    "The [EMIT_L2B_MIN](https://doi.org/10.5067/EMIT/EMITL2BMIN.001) product is generated to support the EMIT mission objectives of constraining the sign of dust related radiative forcing. Ten mineral types are the core focus of this work: Calcite, Chlorite, Dolomite, Goethite, Gypsum, Hematite, Illite+Muscovite, Kaolinite, Montmorillonite, and Vermiculite. A future product will aggregate these results for use in Earth System Models. Additional minerals are included in this product for transparency but were not the focus of this product. Further validation is required to use these additional mineral maps, particularly in the case of resource exploration. Similarly, the separation of minerals with similar spectral features, such as a fine-grained goethite and hematite, is an area of active research. The results presented here are an initial offering, but the precise categorization is likely to evolve over time, and the limits of what can and cannot be separated on the global scale is still being explored. The user is encouraged to read the [Algorithm Theoretical Basis Document (ATBD)](https://lpdaac.usgs.gov/documents/1659/EMITL2B_ATBD_v1.pdf) for more details.\n",
    "\n",
    "**Requirements** \n",
    " - Set up Python Environment - See **setup_instructions.md** in the `/setup/` folder \n",
    "\n",
    "**Learning Objectives**  \n",
    "- How to open an EMIT L2B `.nc` file as an `xarray.Dataset`\n",
    "- Apply the Geometry Lookup Table (GLT) to orthorectify the image.\n",
    "- Find minerals of interest within a granule\n",
    "- Visualize Mineral Identification and Band depth\n",
    "- Evaluate mineral uncertainty\n",
    "- Calculate and Visualize mineral Abundance\n",
    "\n",
    "**Tutorial Outline**  \n",
    "\n",
    "1. Setup  \n",
    "2. Working with Mineral Identification and Band Depth\n",
    "3. Aggregating and Mineral Abundance\n",
    "4. Exporting to Cloud-Optimized GeoTIFF (COG)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b6c3019-765d-47b7-8536-5c24ff0288c4",
   "metadata": {},
   "source": [
    "## 1. Setup\n",
    "\n",
    "Import the required Python libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5370bd60-e09c-4ba2-82f7-09697360ff79",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import earthaccess\n",
    "import geopandas as gp\n",
    "import os\n",
    "import sys\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import hvplot.xarray\n",
    "import holoviews as hv\n",
    "import panel as pn\n",
    "sys.path.append('../modules/')\n",
    "import emit_tools as et"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b3798c2",
   "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,
   "id": "909ddf57",
   "metadata": {},
   "outputs": [],
   "source": [
    "earthaccess.login(persist=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b94d9b18",
   "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 L2B Mineral Identification and Band Depth granule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b858790",
   "metadata": {},
   "outputs": [],
   "source": [
    "# List the browse images from the text file output of the previous notebook.\n",
    "min_list = '../../data/results_urls.txt'\n",
    "with open(min_list) as f:\n",
    "    min_urls = [line.rstrip('\\n') for line in f]\n",
    "min_urls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57a34e59",
   "metadata": {},
   "outputs": [],
   "source": [
    "# List the browse images from the text file output of the previous notebook.\n",
    "unc_list = '../../data/min_uncert_urls.txt'\n",
    "with open(unc_list) as f:\n",
    "    unc_urls = [line.rstrip('\\n') for line in f]\n",
    "unc_urls"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84d63672",
   "metadata": {},
   "source": [
    "Get an HTTPS Session using your earthdata login, set a local path to save the file, and download the granule asset, in this case we are selecting the second scene (index 1) from the list of granules for both the mineral and uncertainty files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb62adee",
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = earthaccess.get_fsspec_https_session()\n",
    "fp = fs.open(min_urls[0])\n",
    "fp_un = fs.open(unc_urls[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4e3769f",
   "metadata": {},
   "source": [
    "## 1.2 Downloaded Data\n",
    "\n",
    "If you’ve already downloaded the data using the workflow shown in Section 6 of the [Finding EMIT L2B Mineralogy Data](Finding_EMIT_L2B_Mineralogy_Data.ipynb) , you can just set filepaths using the cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14114c12-c896-43d9-9424-3d015b98a0b5",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "#fp = '../../data/EMIT_L2B_MIN_001_20230427T173309_2311711_010.nc' # Mineral\n",
    "#fp_rgb = '../../data/EMIT_L2A_RFL_001_20230427T173309_2311711_010.png' # RGB\n",
    "#fp_un = '../../data/EMIT_L2B_MINUNCERT_001_20230427T173309_2311711_010.nc' # Mineral Uncertainty"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c28b8fb",
   "metadata": {},
   "source": [
    "## 2. Working with the L2B Mineral Identification and Band Depth"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39339147",
   "metadata": {},
   "source": [
    "EMITL2BMIN data are distributed in a non-orthocorrected spatially raw NetCDF4 (.nc) format consisting of the data and its associated metadata. Inside the `.nc` file there are 3 groups. Groups can be thought of as containers to organize the data. \n",
    "\n",
    "1. The root group that can be considered the main dataset contains 4 data variables data described by the downtrack, and crosstrack dimensions. These variables are `group_1_mineral_id`, `group_1_band_depth`, `group_2_mineral_id`, and `group_2_band_depth`. These contain the ID and a band depth for each mineral group. These groups do not correspond to the `.netcdf` file groups, but rather the spectral library groups used to identify the minerals based on which region of the spectra the mineral features correspond to.\n",
    "2. The `mineral_metadata`  group containing the spectral library entry name, index, record, group, and url for each entry.\n",
    "3. The `location` group contains latitude and longitude values at the center of each pixel described by the crosstrack and downtrack dimensions, as well as a geometry lookup table (GLT) described by the ortho_x and ortho_y dimensions. The GLT is an orthorectified image (EPSG:4326) consisting of 2 layers containing downtrack and crosstrack indices. These index positions allow us to quickly project the raw data onto this geographic grid.\n",
    "\n",
    "To access the `.nc` file, you can use the `netCDF4`, `xarray` libraries, or fuctions from the `emit_tools.py` library. Here we will use the `emit_xarray` function from this library, which will open and organize the data into an easy to work with `xarray.Dataset` object. We can also pass the `ortho=True` argument to orthorectify the data at this stage, but we will start just examining the data to get a better understanding. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e97803b-d221-45b0-b388-a851b5fee96a",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ds_min = et.emit_xarray(fp)\n",
    "ds_min"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a54c4bd3",
   "metadata": {},
   "source": [
    "If we look at the mineral `index` by printing the first 5 values, we can see that values start with 1. If we look at the minimum values of the mineral IDs we can see these have 0 as a possible value. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6bdd201e",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(ds_min.index.data[:5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "768344b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f'Group_1_minimum:{ds_min.group_1_mineral_id.data.min()} Group_2_minimum: {ds_min.group_2_mineral_id.data.min()}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c98dfb5",
   "metadata": {},
   "source": [
    "The 0 here represents no match.  For convenience, let's make a DataFrame that holds the mineral data, and add that 'No match' reference to it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88e8e8bc-1eb6-4b8c-96ee-2c043b8ce9cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "min_df = pd.DataFrame({x: ds_min[x].values for x in [var for var in ds_min.coords if 'mineral_name' in ds_min[var].dims]})\n",
    "min_df.loc[-1] = {'index': 0, 'mineral_name': 'No_Match', 'record': -1.0, 'url': 'NA', 'group': 1.0, 'library': 'NA', 'spatial_ref': 0}\n",
    "min_df = min_df.sort_index().reset_index(drop=True)\n",
    "min_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e56c228d",
   "metadata": {},
   "source": [
    "### 2.1 Orthorectification\n",
    "\n",
    "The orthorectifation process has already been done for EMIT data. Here we are just using the crosstrack and downtrack indices contained in the GLT to place our spatially raw mineralogy data a into geographic grid with the `ortho_x` and `ortho_y` dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c444096a",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_min = et.ortho_xr(ds_min)\n",
    "ds_min"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe4aa0c5",
   "metadata": {},
   "source": [
    "We can see from these outputs that the dimensions are now latitude and longitude."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a44790aa",
   "metadata": {},
   "source": [
    "In this example, we'll just work with the group_1 mineral data. We can find the minerals present in the scene by finding unique values in the `group_1_mineral_id`, but first we will replace fill-values introduced during orthorectification with `np.nan`, to omit them from our analysis and improve visualizations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d2679d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assign fill to np.nan\n",
    "for var in ds_min.data_vars:\n",
    "    ds_min[var].data[ds_min[var].data == -9999] = np.nan"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3565883f-0d9c-4556-a81f-0adefb971518",
   "metadata": {},
   "source": [
    "## 2.2 Visualize Group 1 Minerals\n",
    "\n",
    "To visualize minerals present, plot Group 1 Minerals using a categorical color set. You can hover over a colored region to see the zero-based mineral id from the spectral library. Note that these values correspond to the 1-based index value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5debcfaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_min['group_1_mineral_id'].hvplot.image(cmap='glasbey', geo=True, tiles='ESRI', alpha=0.8,frame_width=750).opts(title='Group 1 Mineral ID')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fafa07ec",
   "metadata": {},
   "source": [
    "This figure shows the minerals present in the scene, but doesn't really quantify how well they matched with the spectral library. For that we can look at the band depth for each mineral. We can build an interactive tool to do this using the `panel` and `hvplot` libraries. This will take a bit of time to load for each selection.\n",
    "\n",
    "Because many minerals are scarce, we'll start by updating the names to include relative fractions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e490efeb-0e98-4a95-8ed4-9167ee5029c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "g1_min_percent = [np.round(np.sum(ds_min.group_1_mineral_id.data.flatten() == g1min) / np.sum(ds_min.group_1_mineral_id.data.flatten() > 0),2) * 100 for g1min in range(len(min_df))]\n",
    "g1_dropdown_names = [str(g1_min_percent[_x]) + ' %: ' + x for (_x, x) in enumerate(min_df.mineral_name.tolist()) if g1_min_percent[_x] > 0 and x != 'No_Match']\n",
    "g1_dropdown_names = np.array(g1_dropdown_names)[np.argsort([float(x.split(' %:')[0]) for x in g1_dropdown_names])[::-1]].tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b5965751",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Interactive Panel Control For Mineral Band Depth\n",
    "min_select = pn.widgets.Select(name='Mineral Name', options = g1_dropdown_names, value = g1_dropdown_names[0])\n",
    "@pn.depends(min_select)\n",
    "def min_browse(min_select):\n",
    "    mask = ds_min['group_1_band_depth'].where(ds_min['group_1_mineral_id'] == min_df['mineral_name'].tolist().index(min_select.split('%: ')[-1]))\n",
    "    map = mask.hvplot.image(cmap='viridis', geo=True, tiles='ESRI', alpha=0.8,frame_width=450, clim=(0,np.nanpercentile(mask,98))).opts(title=f'{min_select} Band Depth')\n",
    "    return map\n",
    "pn.Row(pn.WidgetBox(min_select),min_browse)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8320cb8-51b2-4f0a-b549-578b0e975dbd",
   "metadata": {},
   "source": [
    "## 2.3 Visualize Group 2 Minerals\n",
    "\n",
    "We can do the same thing with Group 2 Minerals.  Group 2 will show a more diverse set of minerals in this region, including clays and carbonates.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a9862e8-376d-48c8-8af3-6c7a7941326d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds_min['group_2_mineral_id'].hvplot.image(cmap='glasbey', geo=True, tiles='ESRI', alpha=0.8,frame_width=750).opts(title='Group 2 Mineral ID')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8d3110a-5cd8-474c-99c0-70e9e7ca4021",
   "metadata": {},
   "outputs": [],
   "source": [
    "g2_min_percent = [np.round(np.sum(ds_min.group_2_mineral_id.data.flatten() == g1min) / np.sum(ds_min.group_2_mineral_id.data.flatten() > 0),2) * 100 for g1min in range(len(min_df))]\n",
    "g2_dropdown_names = [str(g2_min_percent[_x]) + ' %: ' + x for (_x, x) in enumerate(min_df.mineral_name.tolist()) if g2_min_percent[_x] > 0 and x != 'No_Match']\n",
    "g2_dropdown_names = np.array(g2_dropdown_names)[np.argsort([float(x.split(' %:')[0]) for x in g2_dropdown_names])[::-1]].tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "284282c0-3be6-4a7b-b054-836d3d83ae69",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Interactive Panel Control For Mineral Band Depth\n",
    "min_select_g2 = pn.widgets.Select(name='Mineral Name', options = g2_dropdown_names, value = g2_dropdown_names[0])\n",
    "@pn.depends(min_select_g2)\n",
    "def min_browse_g2(min_select):\n",
    "    # print(min_select)\n",
    "    mask = ds_min['group_2_band_depth'].where(ds_min['group_2_mineral_id'] == min_df['mineral_name'].tolist().index(min_select.split('%: ')[-1]))\n",
    "    map = mask.hvplot.image(cmap='viridis', geo=True, tiles='ESRI', alpha=0.8,frame_width=450, clim=(0,np.nanpercentile(mask,98))).opts(title=f'{min_select} Band Depth')\n",
    "    return map\n",
    "pn.Row(pn.WidgetBox(min_select_g2),min_browse_g2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e5825a2-7865-402f-8e33-cf45a05490fe",
   "metadata": {},
   "source": [
    "## 2.4 Visualizing and filtering with uncertainties\n",
    "\n",
    "Often just as important as visualizing the core mineral detections is displaying the corresponding uncertainties.  The EMIT L2B MIN proudct comes with two different uncertainties (repeated for each group):\n",
    "\n",
    "1. Fit - this is how good of a fit the individual mineral detection was, as estimated by the correlation coefficient of the alignment between the continuum removed library reference and the continuum removed observed spectra (after scaling).\n",
    "2. Band depth uncertainty - this is the propogation of the reflectance uncertainty through the band depth calculation.\n",
    "\n",
    "Let's take a look at each of these by loading the relevant datasets:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a659aac-1e9a-45ae-80e8-42d97f90633b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open an Orthorectify the Uncertainty Data\n",
    "ds_min_unc = et.emit_xarray(fp_un)\n",
    "ds_min_unc = et.ortho_xr(ds_min_unc)\n",
    "ds_min_unc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30c87313-8810-4a29-9358-00e6187d41ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "min_select_fit_g2 = pn.widgets.Select(name='Mineral Name', options = g1_dropdown_names, value = g1_dropdown_names[0])\n",
    "@pn.depends(min_select_fit_g2)\n",
    "def min_browse_fit_g2(min_select):\n",
    "    mask = ds_min_unc['group_1_fit'].where(ds_min['group_1_mineral_id'] == min_df['mineral_name'].tolist().index(min_select.split('%: ')[-1]))\n",
    "    # The fits are scaled down by a factor of two in the current EMIT L2B products, so scale them back up:\n",
    "    mask *= 2\n",
    "    map = mask.hvplot.image(cmap='viridis', geo=True, tiles='ESRI', alpha=0.8,frame_width=450, clim=(0,np.nanpercentile(mask,98))).opts(title=f'{min_select} Fit')\n",
    "    return map\n",
    "pn.Row(pn.WidgetBox(min_select_fit_g2),min_browse_fit_g2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cd9f2fa",
   "metadata": {},
   "source": [
    "## 3. Aggregating and Mineral Abundance\n",
    "\n",
    "The above visualizations walk through the identification of individual library contituents, and visualize band depths.  However, the Tetracorder library used by EMIT contains many substrates that are spectrally distinct, but which may be useful to group together for some science applications. The library also contains many mixtures - both aerial and intimate.\n",
    "\n",
    "To start, the mineral_grouping_matrix from the emit-sds-l2b repository (coppied locally) contains information aggregated from laboratory XRD analyses to attempt to quantify the abundance of different minerals within each constituent.  A -1 in the spreadsheet indicates an unknown but non-zero quantity, which in the few cases in the EMIT-10 columns we assume to be 100%.  Let's open that spreadsheet and take a look:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9eea13cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open Mineral Groupings .csv\n",
    "mineral_groupings = pd.read_csv('../../data/mineral_grouping_matrix_20230503.csv')\n",
    "# The EMIT 10 Minerals are in columns 6 - 17.  Columns after 17 are experimental, and we'll drop for this tutorial:\n",
    "mineral_groupings = mineral_groupings.drop([x for _x, x in enumerate(mineral_groupings) if _x >= 17], axis=1)\n",
    "\n",
    "# Retrieve the EMIT 10 Mineral Names from Columns 7-16 (starting with 0) in .csv\n",
    "mineral_names = [x for _x, x in enumerate(list(mineral_groupings)) if _x > 6 and _x < 17]\n",
    "# Use EMIT 10 Mineral Names to Subset .csv to only columns with EMIT 10 mineral_names\n",
    "mineral_abundance_ref = np.array(mineral_groupings[mineral_names])\n",
    "# Replace Some values in the .csv\n",
    "mineral_abundance_ref[np.isnan(mineral_abundance_ref)] = 0\n",
    "mineral_abundance_ref[mineral_abundance_ref == -1] = 1\n",
    "\n",
    "mineral_groupings"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a560c2ea-26f0-4a6c-9673-3ea65feab1eb",
   "metadata": {},
   "source": [
    "# 3.1 Approximating Mineral Spectral Abundance\n",
    "\n",
    "If we make the assumption that the XRD analysis are accurate, and that band-depth scales linearly with abundance, we can approximate the mineral abundance at the surface.  It should be noted that both of these assumptions are fraught - XRD analyses break down for some very small grainsize particles, particularly the Iron Oxides (Goethite and Hematite), and band-depth is heavily influenced by particle grain size, and so the abundance-band depth relationship is not linear.  We will expand on these details a bit more later, but for now lets take a look at what happens if we hold both of these as true.\n",
    "\n",
    "The first step is to run through each mineral that has a non-nan value in the mineral_groupings dataframe, and add up the sum of each of those within the scene.  A little matrix multiplication is all we need to do that.  Notably, in the emit-sds/emit-sds-l2b, this functionality is already built into the group_aggregator.py script in an efficient manner...but because the calculation is simple we'll reproduce here for learning purposes.\n",
    "\n",
    "Below, we step through each mineral one at a time, and multiply the band depth by the XRD value for each constintuent in the scene, storing it in a (y,x,band) numpy array, which we can cast into x-array for visualization down the line."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b830cda5-ddef-40c3-a1ca-a4dcd310c3d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "mineral_abundance = np.zeros((ds_min['group_1_band_depth'].shape[0], ds_min['group_1_band_depth'].shape[1], len(mineral_names)))\n",
    "for _m, mineral_name in enumerate(mineral_names):\n",
    "    print(f'Calculating {mineral_name}')\n",
    "    for _c in range(mineral_groupings.shape[0]):\n",
    "        if np.isnan(mineral_groupings[mineral_name][_c]) == False:      \n",
    "            group = mineral_groupings[\"Group\"][_c]\n",
    "            mineral_abundance[...,_m] += (ds_min[f'group_{group}_mineral_id'].values == mineral_groupings['Index'][_c]) * ds_min[f'group_{group}_band_depth'].values\n",
    " \n",
    "            \n",
    "mineral_abundance[np.isnan(mineral_abundance)] = 0\n",
    "mineral_abundance[np.all(mineral_abundance == 0, axis=-1),:] = np.nan\n",
    "\n",
    "# Cast as x-array for consistent visualization\n",
    "mineral_abundance_xarray = xr.merge([xr.DataArray(mineral_abundance[...,_x],\n",
    "                                                  name=mineral_names[_x],\n",
    "                                                  coords=ds_min['group_1_band_depth'].coords,\n",
    "                                                  attrs=ds_min.attrs,) \n",
    "                                     for _x in range(len(mineral_names))])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7445a29b-fcae-4c57-a505-0336f0286590",
   "metadata": {},
   "source": [
    "After which we can show them as a widget:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ef3cf63-4fa1-4ac9-a43c-cd023769b8ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "min_abun_wid = pn.widgets.Select(name='Mineral Abundance', options = mineral_names, value = mineral_names[0])\n",
    "@pn.depends(min_abun_wid)\n",
    "def min_abun_browse(mineral_name):\n",
    "    map = mineral_abundance_xarray[mineral_name].hvplot.image(cmap='viridis', geo=True, tiles='ESRI', alpha=0.8,frame_width=450, \n",
    "                                                              clim=(0,np.nanpercentile(mineral_abundance_xarray[mineral_name],98))).opts(\\\n",
    "                                                              title=f'{mineral_name} Spectral Abundance')\n",
    "\n",
    "    return map\n",
    "pn.Row(pn.WidgetBox(min_abun_wid),min_abun_browse)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "326609dd-b759-4e73-b00a-32f082ba2e3e",
   "metadata": {},
   "source": [
    "We can also take a look at the histogram distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27283418-47fb-423e-b848-fda2754d0506",
   "metadata": {},
   "outputs": [],
   "source": [
    "max_bin_max = np.max([np.nanpercentile(mineral_abundance_xarray[x].values,98) for x in mineral_names])\n",
    "mineral_abundance_xarray.hvplot.hist(y=mineral_names, cmap='glasbey', bins=50, alpha=0.5, legend='left', bin_range=(0.001, max_bin_max), \n",
    "                                     xlabel='Spectral Abundance', ylabel='Pixel Count', title='',frame_width=650)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8904ef13-b867-4723-a289-774c318ffa73",
   "metadata": {},
   "source": [
    "Notably, these abundance fractions do not sum to one.  This can either be due to copmlex mixtures that aren't fully captured, or due to the presence of additional materials that do not directly absorb in the VSWIR - most namely quartz and feldspar, but also including organic soil matter, as well as NPV and PV.  To illustrate this, let's examine the totals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc40f552-8830-4e50-b55c-83b994fb5a4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "total_mineral_abundance = xr.DataArray(np.nansum(mineral_abundance,axis=-1), \n",
    "                                       name='Total Abundance', \n",
    "                                       coords=ds_min['group_1_band_depth'].coords, \n",
    "                                       attrs=ds_min.attrs,)\n",
    "total_mineral_abundance.hvplot.hist(y=['Total Abundance'], cmap='glasbey', bins=50, alpha=0.5, legend='left', \n",
    "                                    bin_range=(0.001, np.nanpercentile(total_mineral_abundance.values, 99.9)), \n",
    "                                    xlabel='Spectral Abundance', ylabel='Pixel Count', title='',frame_width=650)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0fb01dad",
   "metadata": {},
   "source": [
    "## 4. Exporting To Cloud-Optimized Geotiffs\n",
    "\n",
    "We can select layers from any of the xarray datasets we've created and export them to a Cloud-Optimized Geotiff (COG) using the `rio.to_raster` function from the `rasterio` library. This will allow us to share the data with others or use it in other GIS software.\n",
    "\n",
    "First create an output folder."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "677f044d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create Out filenames and set folder\n",
    "out_folder = '../../data/output/' # may need to change based on your directory structure\n",
    "# Create out_folder if it does not exist\n",
    "if not os.path.exists(out_folder):\n",
    "   os.makedirs(out_folder)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e201ddd",
   "metadata": {},
   "source": [
    "### 4.1 Mineral ID"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a83fd5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set output Filename\n",
    "out_name = f'{ds_min.granule_id}_group_1_mineral_id.tif'\n",
    "# Select Group to Output\n",
    "dat_out = ds_min['group_1_mineral_id']\n",
    "# Replace nan with a fill value\n",
    "dat_out.data = np.nan_to_num(dat_out.data, nan=-9999)\n",
    "# Change datatype for integers/categorical\n",
    "dat_out.data = dat_out.data.astype(int)\n",
    "# Encode nodata value\n",
    "dat_out.rio.write_nodata(-9999, encoded=True, inplace=True)\n",
    "# Write data to COG\n",
    "dat_out.rio.to_raster(raster_path=f'{out_folder}{out_name}', driver='COG')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8499033",
   "metadata": {},
   "source": [
    "### 4.2 Band Depth"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f32807f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set output Filename\n",
    "out_name = f'{ds_min.granule_id}_group_1_band_depth.tif'\n",
    "# Select Group to Output\n",
    "dat_out = ds_min['group_1_band_depth']\n",
    "# Replace nan with a fill value\n",
    "dat_out.data = np.nan_to_num(dat_out.data, nan=-9999)\n",
    "# Encode nodata value\n",
    "dat_out.rio.write_nodata(-9999, encoded=True, inplace=True)\n",
    "# Write data to COG\n",
    "dat_out.rio.to_raster(raster_path=f'{out_folder}{out_name}', driver='COG')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f3cbdb50",
   "metadata": {},
   "source": [
    "### 4.3 Abundance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f62e0222",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set output Filename\n",
    "out_name = f'{mineral_abundance_xarray.granule_id}_Calcite.tif'\n",
    "# Select Layer\n",
    "dat_out = mineral_abundance_xarray['Calcite']\n",
    "# Replace nan with a fill value\n",
    "dat_out.data = np.nan_to_num(dat_out.data, nan=-9999)\n",
    "# Encode nodata value\n",
    "dat_out.rio.write_nodata(-9999, encoded=True, inplace=True)\n",
    "# Write data to COG\n",
    "dat_out.rio.to_raster(raster_path=f'{out_folder}{out_name}', driver='COG')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9996cfae",
   "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: 07-07-2024  \n",
    "\n",
    "¹Work performed under USGS contract 140G0121D0001 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.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}