Mineral Indices

Background: Landsat 8 data is available from April 2013 onwards. Mineral Indices can be used to aid in geological mapping, and exploit differences in the spectral response of landscape and regolith materials. Care needs to be taken in the interpretation of these indices as the scale of the imagery and sensor characteristics can affect the usefulness of the indices.

Ratios in this notebook:

\begin{align} Iron\ Oxide\ Ratio\ &= \frac{RED}{BLUE} \\ \end{align}
\begin{align} Ferrous\ Minerals\ Ratio &= \frac{SWIR1}{NIR} \\ \end{align}
\begin{align} Clay\ NIR/SWIR1\ Ratio &= \frac{NIR}{SWIR1} \\ \end{align}
\begin{align} Clay\ Minerals\ Ratio &= \frac{SWIR1}{SWIR2} \\ \end{align}
\begin{align} Ferruginous\ Regolith\ Ratio &= \frac{NIR}{GREEN} \\ \end{align}

Equation (1) Iron Oxide Ratio: (Drury 1987, Harris Geospatial Solutions,Inc. 2018, Segal 1982)

Equation (2) Ferrous Minerals Ratio: (Drury 1987, Harris Geospatial Solutions,Inc. 2018, Segal 1982)

Equation (3) Clay NIR/SWIR Ratio: Reference not found

Equation (4) Clay Minerals Ratio: (Drury 1987, Harris Geospatial Solutions,Inc. 2018, Wilford and Creasey 2002)

Equation (5) Ferruginous Regolith Ratio: (Wilford and Creasey 2002)

Before you run this notebook: You need to run the following commands from the command line prior to launching jupyter notebook from the same terminal so that the required libraries and paths are set.

module use /g/data/v10/public/modules/modulefiles

module load dea

What does this notebook do?: This notebook runs mineral indices on Landsat 8 data. You can then save the data as a png image, netcdf and geotiff.

Date: August 2018

Authors: Leo Lymburner, Bex Dunn

Tags: Landsat8, Mineral IndicesGeotiff, NetCDF

Import modules from standard libraries, datacube and files

Select ‘Trust this notebook’ to import these modules.

[1]:
%pylab notebook

import os

# modules for datacube
import datacube
from datacube.helpers import write_geotiff

# set datacube alias (just a string with what you're doing)
dc = datacube.Datacube(app='dc-Mineral Indices')

# Import external functions from dea-notebooks
sys.path.append(os.path.expanduser('~/dea-notebooks/10_Scripts/'))
import BandIndices
import DEADataHandling
Populating the interactive namespace from numpy and matplotlib

Edit this cell to choose your area and temporal range

[2]:
#Define temporal range
start_of_epoch = '2013-01-01'
end_of_epoch =  '2013-06-12'

#Use this to manually define an upper left/lower right coords
lat_max = -15.2
lat_min = -15.5
lon_max = 128.3
lon_min = 128.6

query = {'time': (start_of_epoch, end_of_epoch),
         'x' : (lon_min, lon_max),
         'y' : (lat_max, lat_min),
         'crs' : 'EPSG:4326'} #using WGS 84 here

Load data from Digital Earth Australia

[3]:
#call our retrieved xarray dataset ls8_nbart (because we're using nbart see DEA link in heading above)
ls8_nbart, crs, affine = DEADataHandling.load_nbarx(dc=dc, sensor='ls8',query=query)
Loading ls8_nbart_albers
Loaded ls8_nbart_albers
Generating mask ls8_pq_albers
[4]:
ls8_nbart
[4]:
<xarray.Dataset>
Dimensions:          (time: 3, x: 1336, y: 1355)
Coordinates:
  * time             (time) datetime64[ns] 2013-04-26T01:31:31 ...
  * y                (y) float64 -1.618e+06 -1.618e+06 -1.618e+06 -1.618e+06 ...
  * x                (x) float64 -4.007e+05 -4.007e+05 -4.007e+05 -4.007e+05 ...
Data variables:
    coastal_aerosol  (time, y, x) float64 544.0 581.0 560.0 524.0 520.0 ...
    blue             (time, y, x) float64 508.0 535.0 511.0 481.0 474.0 ...
    green            (time, y, x) float64 708.0 730.0 700.0 652.0 665.0 ...
    red              (time, y, x) float64 919.0 952.0 907.0 862.0 870.0 ...
    nir              (time, y, x) float64 1.279e+03 1.299e+03 1.264e+03 ...
    swir1            (time, y, x) float64 821.0 900.0 891.0 794.0 821.0 ...
    swir2            (time, y, x) float64 373.0 464.0 445.0 367.0 367.0 ...
Attributes:
    crs:      EPSG:3577
    affine:   | 25.00, 0.00,-400750.00|\n| 0.00,-25.00,-1618400.00|\n| 0.00, ...

Calculate Mineral Indices

\begin{align} Iron\ Oxide\ Ratio\ &= \frac{RED}{BLUE} \\ \end{align}
\begin{align} Ferrous\ Minerals\ Ratio &= \frac{SWIR1}{NIR} \\ \end{align}
\begin{align} Clay\ NIR/SWIR1\ Ratio &= \frac{NIR}{SWIR1} \\ \end{align}
\begin{align} Clay\ Minerals\ Ratio &= \frac{SWIR1}{SWIR2} \\ \end{align}
\begin{align} Ferruginous\ Regolith\ Ratio &= \frac{NIR}{GREEN} \\ \end{align}
[5]:
iron_oxide_ratio =  BandIndices.geological_indices(ls8_nbart, 'IOR') #Iron Oxide Ratio (RED/BLUE)
ferrous_iron =  BandIndices.geological_indices(ls8_nbart, 'FMR') #Ferrous Minerals Ratio (SWIR1/NIR)
clay_nirswir = ls8_nbart.nir/ls8_nbart.swir1  #Clay NIR/SWIR1 (NIR/SWIR1)
clay_swir1_swir2 = BandIndices.geological_indices(ls8_nbart, 'CMR') #Clay Minerals Ratio (SWIR1/SWIR2)
ferruginous_reg = ls8_nbart.nir/ls8_nbart.green #Ferruginous Regolith Ratio (NIR/GREEN)
The formula we are using for Iron Oxide Ratio is (red / blue)
The formula we are using for Ferrous Minerals Ratio is (swir1 / nir)
The formula we are using for Clay Minerals Ratio is (swir1 / swir2)

Plot the Iron Oxide Ratio (RED/BLUE)

[6]:
iron_oxide_ratio
[6]:
<xarray.DataArray (time: 3, y: 1355, x: 1336)>
array([[[1.809055, 1.779439, ...,      nan,      nan],
        [1.724395, 1.725455, ...,      nan,      nan],
        ...,
        [1.333966, 1.431686, ..., 1.563969, 1.508242],
        [1.363806, 1.423019, ..., 1.5     , 1.526455]],

       [[1.807829, 1.831239, ..., 1.802057, 1.81491 ],
        [1.823105, 1.830909, ..., 1.78628 , 1.790404],
        ...,
        [1.397647, 1.483165, ..., 1.661836, 1.634085],
        [1.322506, 1.472628, ..., 1.692118, 1.650873]],

       [[     nan,      nan, ..., 1.714894, 1.710638],
        [     nan,      nan, ..., 1.683544, 1.705637],
        ...,
        [1.324385, 1.388031, ..., 1.54321 , 1.49359 ],
        [1.283063, 1.423459, ..., 1.563786, 1.52795 ]]])
Coordinates:
  * time     (time) datetime64[ns] 2013-04-26T01:31:31 ...
  * y        (y) float64 -1.618e+06 -1.618e+06 -1.618e+06 -1.618e+06 ...
  * x        (x) float64 -4.007e+05 -4.007e+05 -4.007e+05 -4.007e+05 ...
[7]:
iron_oxide_ratio_stats = iron_oxide_ratio.max(dim = 'time') #you can change .median to .min, .mean, .max
iron_oxide_ratio_stats.plot( vmin = 0, vmax = 3)
plt.title('Maximum Iron Oxide Ratio (RED/BLUE)')
plt.show()

Plot the Ferrous Minerals Ratio (SWIR1/NIR)

[8]:
ferrous_iron_stats = ferrous_iron.max(dim = 'time') #you can change .median to .min, .mean, .max
fig = plt.figure()
ferrous_iron_stats.plot(vmin = 0, vmax  = 2)
plt.title('Maximum Ferrous Minerals Ratio (SWIR1/NIR)')
plt.show()