AVIRIS
AVIRIS Hyperspectral Sensor
Introduction
AVIRIS is an acronym for the Airborne Visible InfraRed Imaging Spectrometer (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 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:
gain = []
with open(AVIRISDir+"rdn_e/"+AVIRISDir+"rdn_e.gain") as file:
for line in file:
gain.append(line.strip().split()[0])
After that, the DN to Radiance correction can be done like this:
# 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()
Import in GRASS, gain correction (DN2Rad) & atmospheric correction (Rad2Ref)
#!/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