View source for AVIRIS

== AVIRIS Hyperspectral Sensor ==
[[File:aviris_text2.png]]
=== Introduction ===
'''AVIRIS''' is an acronym for the Airborne Visible InfraRed Imaging Spectrometer ([https://aviris.jpl.nasa.gov/ AVIRIS Main page]). 

Its optical sensor creates calibrated images of the upwelling spectral radiance in 224 contiguous spectral channels/bands with wavelengths from 400 to 2500 nanometers (nm). 

=== Airplanes ===
'''AVIRIS''' has been flown on four aircraft platforms: 
* NASA's ER-2 jet
* Twin Otter International's turboprop 
* Scaled Composites' Proteus 
* NASA's WB-57. 

Flight details:
* The ER-2 flies at approximately 20km above sea level, at about 730 km/hr. 
* The Twin Otter aircraft flies at 4km above ground level at 130km/hr. 

=== Data Download ===
Downloading data from '''AVIRIS''' flights is straightforward using the [https://aviris.jpl.nasa.gov/dataportal/ AVIRIS Data Portal].

== Import AVIRIS in GRASS ==
=== Data Downloaded ===

'''AVIRIS''' is downloaded as a .tar.gz file (i.e. ''f180116t01p00r07.tar.gz'') uncompressing as a directory (i.e. ''f180116t01p00r07rdn_e/'') containing with a set files like this:

==== List of files ====
* AVIRIS_OrthoProcessing_Info.txt
* f180116t01p00r07rdn_e_eph
* f180116t01p00r07rdn_e.gain
* f180116t01p00r07rdn_e_lonlat_eph
* ''f180116t01p00r07rdn_e_obs''
* f180116t01p00r07rdn_e_obs.hdr
* ''f180116t01p00r07rdn_e_obs_ort''
* f180116t01p00r07rdn_e_obs_ort.hdr
* ''f180116t01p00r07rdn_e_ort_glt''
* f180116t01p00r07rdn_e_ort_glt.hdr
* ''f180116t01p00r07rdn_e_ort_igm''
* f180116t01p00r07rdn_e_ort_igm.hdr
* f180116t01p00r07rdn_e_ort.plog
* f180116t01p00r07rdn_e.rccf18
* ''f180116t01p00r07rdn_e_sc01_ort_img''
* f180116t01p00r07rdn_e_sc01_ort_img.hdr
* f180116t01p00r07rdn_e.spc

All files twinned with a .hdr file, and in italics in the list above, are ENVI style binary imagery datasets (''BSQ format''). Towards the end of the list, the file matching the pattern ''sc01'' is the largest file, and the one holding the actual hyperspectral data.

==== DN to Radiance Conversion ====

The DN to Radiance correction is done by applying a gain to the stored data in the BSQ file. In this case, the gain file in ''f180116t01p00r07rdn_e.gain''. Extracting the gain values can be done like this in Python:

<syntaxhighlight lang="python">
gain = []
with open(AVIRISDir+"rdn_e/"+AVIRISDir+"rdn_e.gain") as file:
        for line in file:
            gain.append(line.strip().split()[0])
</syntaxhighlight>

After that, the DN to Radiance correction can be done like this:

<syntaxhighlight lang="python">
# define GRASS-Python environment
os.environ['GISBASE'] = "/usr/local/grass84"
grass_pydir = os.path.join("/usr/local/grass84", "etc", "python")
sys.path.append(grass_pydir)
from grass.script import setup as gsetup
from grass.pygrass.modules.shortcuts import raster as r
from grass.pygrass.modules import Module

rmapcalc = Module('r.mapcalc')
# AVIRIS DN->Rad->Ref
    print("Convert DN to Rad")
    # Apply Gain from file
    for j in range(len(gain)):
        exp = str("rad."+str(j+1)+"=dn."+str(j+1)+"*10.0/"+str(gain[j]))
        #print(exp)
        rmapcalc(expression=exp, overwrite=OVR, quiet=QUIET)
        rmapcalc.wait()
</syntaxhighlight>

=== Import in GRASS, gain correction (DN2Rad) & atmospheric correction (Rad2Ref) ===
<syntaxhighlight lang="python">
#!/usr/bin/env python
"""
Process all AVIRIS radiance tar.gz into reflectance.

Download all radiance tar.gz files into rsdatapath
Make sure GRASS GIS and python-gdal are installed
run this code after specifying your directories

# Date: 24th October, 2023
# Public domain, GRASS GIS

# PURPOSE
# This script processes AVIRIS images
# 1 - untar *.tar.gz files
# 2 - Connect files in GRASS GIS Location of your choice
# 3 - DN to Radiance (Apply Gain)
# 4 - Radiance to Surface reflectance (i.atcorr)
"""
# Load necessary libraries
import os
import sys
import time
import struct
from subprocess import PIPE
# from osgeo import ogr
from osgeo import osr
from osgeo import gdal
# define GRASS-Python environment
os.environ['GISBASE'] = "/usr/local/grass84"
grass_pydir = os.path.join("/usr/local/grass84", "etc", "python")
sys.path.append(grass_pydir)
from grass.script import setup as gsetup
from grass.pygrass.modules.shortcuts import raster as r
from grass.pygrass.modules.shortcuts import imagery as i
from grass.pygrass.modules import Module
from grass.script.utils import parse_key_val

# QUIET REPORTS
QUIET = True

# OVERWRITE EXISTING FILES
OVR = True

# Setup the path to the AVIRIS Directories
rsdatapath = "/home/yann/RSDATA/AVIRIS"

# Setup your GRASS GIS working directory
grassdata = "/home/yann/grassdata/"
location = "AVIRIS"
print(location)
mapset = "PERMANENT"

# DEM input for atmospheric correction
inDEM = "/home/yann/RSDATA/SRTM/dem.vrt"

# Visibility distance [Km]
vis = 18

# List of .tar.gz AVIRIS radiance tracks downloaded
AVIRISDirs = ["f180116t01p00r07"]  # , "f100517t01p00r11"]

# CONFIGURATION OF 6S
# https://grass.osgeo.org/grass84/manuals/i.atcorr.html

# Atmospheric mode
# none = 0
# tropical = 1
# midlatitude summer = 2
# midlatitude winter = 3
atm_mode = 1
# Aerosol model
# none = 0
# continental = 1
# maritime = 2
# urban = 3
aerosol_mode = 2

# END OF USER CHANGES
# DO NOT CHANGE ANYTHING AFTER THIS LINE !

# Change directory to find the tar.gz files of AVIRIS
os.chdir(rsdatapath)

# Remove previous temporary GRASS Location
p = os.system("rm -Rf "+grassdata+location)


def read_ephemeris_file(file_path):
    """
    Read the ephemeris _eph file and extract flight altitude.

    Parameters
    ----------
    file_path : STRING
        the ephemeris *_eph file.

    Returns
    -------
    data : FLOAT
        the average height of sensor during flight.
    """
    data = []
    # Define the format of a single record in the binary file
    record_format = 'd' * 6  # Six double-precision floating-point numbers
    with open(file_path, 'rb') as f:
        while True:
            record = f.read(struct.calcsize(record_format))
            if not record:
                break  # Reached end of file
            values = struct.unpack(record_format, record)
            data.append(values)
    return data


# START PROCESS
# Read Gain from File
gain = []
# Start processing each radiance track bundlefrom DN to Ref
for AVIRISDir in AVIRISDirs:
    # Untar all of your AVIRIS images in all your directories
    print("Untar AVIRIS files in ", AVIRISDir)
    #p = os.system("tar xvf "+AVIRISDir+".tar.gz")
    with open(AVIRISDir+"rdn_e/"+AVIRISDir+"rdn_e.gain") as file:
        for line in file:
            gain.append(line.strip().split()[0])
    # Process a North alignment of the original file TODO add datatype=float
    kwargs = {'format': 'GTiff', 'geoloc': True, 'options': 'overwrite'}
    # gdal.UseExceptions()
    #ds = gdal.Warp("del.tif", AVIRISDir + "rdn_e/" + AVIRISDir + "rdn_e_sc01_ort_img", **kwargs)
    #del ds
    # Create a temporary GRASS Location
    p = os.system("grass -c del.tif -e "+grassdata+location)
    time.sleep(20)
    session = gsetup.init(grassdata, location, mapset)
    print("Import in GRASS GIS")
    # Import DEM for atmospheric correction
    ringdal = Module('r.in.gdal')
    # In case the DEM file is not empty
    # ringdal(input=AVIRISDir+"rdn_e/"+AVIRISDir+"rdn_e_ort_igm", output="dem",
    #         overwrite=OVR)
    # In case the DEM is available to import from srtm.csi.cgiar.org
    #rimport = Module('r.import')
    #rimport(input=inDEM, output="dem", extent="region", overwrite=OVR)
    # The least effective, make your own flat DEM layer
    rmapcalc = Module('r.mapcalc')
    rmapcalc(expression="dem=1", overwrite=OVR)
    rmapcalc.wait()
    # Create a visibility (AOD) map
    rmapcalc(expression="vis=18", overwrite=OVR)
    rmapcalc.wait()
    # Import AVIRIS in GRASS
    ringdal(input="del.tif", output="dn", quiet=QUIET)
    # Remove nodata from DN maps
    rnull = Module('r.null')
    for idx in range(224):
        rnull(map="dn."+str(idx+1), setnull=[-50, 0], quiet=QUIET)
    # Read rotated AVIRIS image
    src = gdal.Open("del.tif")
    ulx_, xres, xskew, uly_, yskew, yres = src.GetGeoTransform()
    lrx_ = ulx_ + (src.RasterXSize * xres)
    lry_ = uly_ + (src.RasterYSize * yres)
    # Setup the source projection
    source = osr.SpatialReference()
    source.ImportFromWkt(src.GetProjection())
    del src
    # The target projection is ll wgs84
    target = osr.SpatialReference()
    target.ImportFromEPSG(4326)
    # Create the transform
    transform = osr.CoordinateTransformation(source, target)
    # Transform the points
    uly, ulx, _ = transform.TransformPoint(ulx_, uly_)
    lry, lrx, _ = transform.TransformPoint(lrx_, lry_)
    # Import Auxiliary layers, rotate it first
    ds = gdal.Warp("del1.tif", AVIRISDir + "rdn_e/" + AVIRISDir +
                   "rdn_e_obs_ort", **kwargs)
    del ds
    # Import in GRASS TODO force Float data
    ringdal(input="del1.tif", output="obs", flags="e", overwrite=OVR, quiet=QUIET)
    # Remove nodata from UTC map
    rnull(map="obs.10", setnull=[-9999, 0], quiet=QUIET)
    # Extract stats
    rinfo = Module('r.info', map="obs.10", flags="s", stdout_=PIPE)
    rinfo_out = parse_key_val(rinfo.outputs.stdout)
    # Extract mean values from Dict
    dectime = rinfo_out["mean"]
    print("timestamp ", str(dectime))
    print(AVIRISDir)
    date = ["20"+AVIRISDir[1:3], AVIRISDir[3:5], AVIRISDir[5:7]]
    print(date)
    # AVIRIS DN->Rad->Ref
    print("Convert DN to Rad")
    # Apply Gain from file
    for j in range(len(gain)):
        exp = str("rad."+str(j+1)+"=dn."+str(j+1)+"*10.0/"+str(gain[j]))
        #print(exp)
        rmapcalc(expression=exp, overwrite=OVR, quiet=QUIET)
        rmapcalc.wait()
    # Atmospheric Correction
    print("Atmospheric Correction")
    # Basic script for i.atcorr for AVIRIS
    # Geometrical conditions (AVIRIS = 31)
    geom = 31
    # Sensor height (satellite is -1000)
    eph_data = read_ephemeris_file(AVIRISDir+"rdn_e/"+AVIRISDir+"rdn_e_eph")
    # Extract average altitude
    add = 0
    for j in range(len(eph_data)):
        add += eph_data[j][2]
    height = add/len(eph_data)
    # Define the sensor height as Km and negative (for 6S)
    sens_height = -1*height/1000.0
    # Here we suppose you have altitude (DEM) and Visibility (VIS) maps ready
    # ---------------------------------------------
    # Altitude dummy value (in Km should be negative in this param file)
    # (overwritten by DEM raster input)
    alt = -0.600
    # datetime of satellite overpass (month, day, GMT decimal hour)
    mdh = str(date[1])+" "+str(date[2])+" "+str("%.3f" % float(dectime))
    print(mdh)
    # Central Lat/Long
    Lon = ulx+(lrx-ulx)/2.0
    Lat = lry+(uly-lry)/2.0
    print(Lon, Lat)
    # satellite band number (AVIRIS [209-432])
    satbandno = 209  # Band 1 of AVIRIS is first for atmospheric correction
    for idx in range(224):
        b = "rad." + str(idx+1)
        b_out = b.replace("rad.", "ref.")
        # Extract min max radiance stats
        rinfo = Module('r.info', map="rad."+str(idx+1), flags="s", stdout_=PIPE)
        rinfo_out = parse_key_val(rinfo.outputs.stdout)
        # Extract mean values from Dict
        minimum = rinfo_out["min"]
        maximum = rinfo_out["max"]
        if (minimum < 0.0):
            minimum = 0.0
        param = []
        # geometrical conditions=AVIRIS
        param.append(str(geom)+"\n")
        # month day hh.ddd longitude latitude (hh.ddd is in decimal hours GMT)
        param.append(str(mdh)+" "+str("%.2f" % Lon)+" "+str("%.2f" % Lat)+"\n")
        # atmospheric mode=tropical
        param.append(str(atm_mode)+"\n")
        # aerosols model=maritime
        param.append(str(aerosol_mode)+"\n")
        # visibility [km] (aerosol model concentration), not used (raster input)\n")
        param.append(str(vis)+"\n")
        # mean target elevation above sea level [km] (here 600m asl), not used (raster input)
        param.append(str(alt)+"\n")
        # sensor height (here, sensor on board a plane)
        param.append(str("%.2f" % round(sens_height,2))+"\n")
        # puw (water vapor content between the aircraft and the surface)
        param.append("-1.0\n")
        # po3 (ozone content between the aircraft and the surface)
        param.append("-1.0\n")
        # taerp (the aerosol optical thickness at 550nm between the aircraft and the surface)
        param.append("-1.0\n")
        # i th band of AVIRIS (atcorr internal no)
        param.append(str(satbandno))
        prm = "param_AVIRIS.txt"
        f = open(prm, "w")
        f.writelines(param)
        f.close()
        print(b)
        print(b_out)
        iatcorr = Module('i.atcorr')
        iatcorr(input=b, elevation="dem", visibility="vis", parameters=prm,
                output=b_out, range=[minimum, maximum], rescale=[0.0,1.0], quiet=QUIET, overwrite=OVR)
        routgdal = Module('r.out.gdal')
        routgdal(input=b_out, output=AVIRISDir+"_"+b_out+".tif", overwrite=OVR, quiet=QUIET)
        satbandno = satbandno + 1
</syntaxhighlight>