Hyperion

From GRASS-Wiki
Jump to navigation Jump to search

Hyperion Hyperspectral Data

EO-1 platform

The NASA EO-1 satellite was launched on November 21, 2000 as part of a one-year technology validation/demonstration mission (EO-1 Mission Page).

The mission was to use the Advanced Land Imager (ALI) instrument on EO-1 to validate and demonstrate technology for the Landsat Data Continuity Mission (LDCM).

The original EO-1 Mission was successfully completed in November 2001. As the end of the Mission approached, the remote sensing research and scientific communities expressed high interest in continued acquisition of image data from EO-1. Based on this user interest and willingness to assist in funding continued operations, an agreement was reached between NASA and the USGS to allow continuation of the EO-1 Program as an Extended Mission.

On February 22nd, 2017, the Earth-Observing One (EO-1) satellite was decommissioned.

Sensors

  • Advanced Land Imager (ALI)
  • Hyperion
  • Linear Etalon Imaging Spectrometer Array (LEISA) Atmospheric Corrector (LAC)

The EO-1 Extended Mission was chartered to collect and distribute ALI multispectral and Hyperion hyperspectral products in response to Data Acquisition Requests (DARs). Under the Extended Mission provisions, image data acquired by EO-1 were archived and distributed by the USGS Center for Earth Resources Observation and Science (EROS) and placed in the public domain.

Hyperion Data Download

Level 1 GST is terrain corrected and provided in 16-bit radiance values. The data are available in Geographic Tagged Image-File Format (GeoTIFF) and are distributed via download at no charge through either EarthExplorer or USGS Global Visualization Viewer (GloVis).

Level 1 GST specifications
Map projection UTM (default zone of the scene center coordinates)
Horizontal Datum WGS84
Resampling method CC (cubic convolution)
Image orientation Map (north up)
Pixel size 30 meter (10 meter pan band)
Format GeoTIFF
Output media Download (no charge)

Extracting and preparing data

The data is downloaded as a zip file (i.e. EO1H1260582015073110T1_1GST.zip) uncompressing as a directory (i.e. EO1H1260582015073110T1_1GST/). Inside exists another directory (i.e. EO1H1260582015073110T1/).

The files inside are following this type of listing:

  • EO1H1260582015073110T1_B001_L1GST.TIF
  • EO1H1260582015073110T1_B002_L1GST.TIF
  • ...
  • EO1H1260582015073110T1_B241_L1GST.TIF
  • EO1H1260582015073110T1_B242_L1GST.TIF
  • EO1H1260582015073110T1_MTL_L1GST.TXT
  • EO1H1260582015073110T1_PF2_01.fgdc
  • README.txt

Of which EO1H1260582015073110T1_MTL_L1GST.TXT is a Landsat type metadata header file (MTL type) well known to Landsat users.

Data perculiarity

It seems that bands 1-8 are filled with NODATA, the same for bands 58 to 76, and again for 225 to 242.

For this reason, the pre-processing script below is avoiding the atmospheric correction of those bands.

This seems to be explained by the presence of two sensors inside the instrument, one dedicated to VNIR hyperspectral sensing and the second one dedicated to SWIR hyperspectral sensing. This is further reinforced by observing the overlapping of band wavelengths across bands 56-57 and 77-78.

Import in GRASS, DN2Rad & Rad2Ref

#!/usr/bin/env python
"""
Process all HYPERION radiance .zip into reflectance.

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

# Date: 1st December, 2023
# Public domain, GRASS GIS

# PURPOSE
# This script processes HYPERION images
# 1 - unzip *.zip 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 glob
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 HYPERION Directories
rsdatapath = "/home/yann/RSDATA/HYPERION"

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

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

# Visibility distance [Km]
vis = 18

# List of .zip HYPERION radiance tracks downloaded
HYPERIONDirs = ["EO1H1260582015073110T1_1GST"]  # , ""]

# 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 HYPERION
os.chdir(rsdatapath)

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

# Define a KVP retrieval function
def findkvp(text, key):
    for item in text.split("\n"):
        if key in item:
            return (item.strip())

# START PROCESS
# Read Gain from File
gain = []
# Start processing each radiance track bundlefrom DN to Ref
for HYPERIONDir in HYPERIONDirs:
    # Unzip all of your HYPERION images in all your directories
    #print("Unzip HYPERION files in ", HYPERIONDir)
    #p = os.system("unzip "+HYPERIONDir+".zip")
    # Go into the new created directory
    path = rsdatapath+'/'+HYPERIONDir+'/'+HYPERIONDir.split('_')[0]
    os.chdir(path)
    # Select the first file in the directory
    first = os.listdir(path)[0]
    print(first)
    # Create a temporary GRASS Location
    p = os.system("grass -c "+first+" -e "+grassdata+location)
    time.sleep(5)
    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 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 (unit is Km)
    rmapcalc = Module('r.mapcalc')
    rmapcalc(expression="dem=0.1", overwrite=OVR)
    rmapcalc.wait()
    # Create a visibility (AOD) map (unit is km)
    rmapcalc(expression="vis=18", overwrite=OVR)
    rmapcalc.wait()
    # Import HYPERION in GRASS
    count = 1
    # To remove nodata from DN maps
    rnull = Module('r.null')
    # Define TIF list
    ext=['TIF','tif']
    listoftiffiles = [f for f in os.listdir(path) if any(f.endswith(e) for e in ext)]
    # Process import
    for file in listoftiffiles:
        outfile='dn.'+str(count)
        ringdal(input=file, output=outfile, quiet=QUIET)
        rnull(map=outfile, setnull=0, quiet=QUIET)
        count = count + 1
    # Read HYPERION image CRS TODO check that "first" has full path attached
    src = gdal.Open(first)
    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_)
    # Read MTL file to import scene characteristics
    MTLfile = glob.glob('*MTL*', recursive=False)
    print("Metadatafile is ",MTLfile)
    with open(MTLfile[0], 'r') as file:
        data = file.read()

    gain_vnir = findkvp(data,"SCALING_FACTOR_VNIR").split("=")[1]
    gain_swir = findkvp(data,"SCALING_FACTOR_SWIR").split("=")[1]
    print("Gains VNIR %d & SWIR %d" % (int(gain_vnir), int(gain_swir)))
    # Extract time mean values
    starttimestamp = findkvp(data,"START_TIME").split("=")[1].strip()
    endtimestamp = findkvp(data,"END_TIME").split("=")[1].strip()
    starttime = starttimestamp.split(" ")[-1]
    endtime = endtimestamp.split(" ")[-1]
    decstarttime = float(starttime.split(":")[0])+float(starttime.split(":")[1])/60.0+float(starttime.split(":")[2])/3600.0
    decendtime = float(endtime.split(":")[0])+float(endtime.split(":")[1])/60.0+float(endtime.split(":")[2])/3600.0
    dectime = 0.5 * (decstarttime + decendtime)
    print("timestamp ", str(dectime))
    print(HYPERIONDir)
    acqdate = findkvp(data, "ACQUISITION_DATE").split("=")[1].strip()
    date = [acqdate.split("-")[0], acqdate.split("-")[1], acqdate.split("-")[2]]
    print(date)
    # HYPERION DN->Rad->Ref
    print("Convert DN to Rad")
    # Apply Gain for VNIR
    for j in range(8, 57+1):
        exp = str("rad."+str(j)+"=dn."+str(j)+"*10.0/"+str(gain_vnir))
        #print(exp)
        rmapcalc(expression=exp, overwrite=OVR, quiet=QUIET)
        rmapcalc.wait()
    # Apply Gain for VNIR
    for j in range(77, 224+1):
        exp = str("rad."+str(j)+"=dn."+str(j)+"*10.0/"+str(gain_swir))
        #print(exp)
        rmapcalc(expression=exp, overwrite=OVR, quiet=QUIET)
        rmapcalc.wait()
    # Atmospheric Correction
    print("Atmospheric Correction")
    # Basic script for i.atcorr for HYPERION
    # Sensor height (satellite is -1000)
    # Define the sensor height as Km and negative (for 6S)
    sens_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)
    # Create the handle for r.out.gdal
    routgdal = Module('r.out.gdal')
    # Process VNIR
    # satellite band number (HYPERION VNIR [433-481])
    satbandno = 433  # Band 8 of HYPERION is first for atmospheric correction
    for idx in range(8,57+1):
        # Geometrical conditions (HYPERION VNIR = 32, SWIR = 33)
        geom = 32
        b = "rad." + str(idx)
        b_out = b.replace("rad.", "ref.")
        # Extract min max radiance stats
        rinfo = Module('r.info', map="rad."+str(idx), flags="s", stdout_=PIPE)
        rinfo_out = parse_key_val(rinfo.outputs.stdout)
        # Extract mean values from Dict
        minimum = float(rinfo_out["min"])
        maximum = float(rinfo_out["max"])
        if (minimum < 0.0):
            minimum = 0.0
        param = []
        # geometrical conditions=HYPERION
        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, a satellite)
        param.append(str("%.2f" % round(sens_height,2))+"\n")
        # i th band of HYPERION (atcorr internal no)
        param.append(str(satbandno))
        prm = "param_HYPERION.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(input=b_out, output=HYPERIONDir+"_"+b_out+".tif", overwrite=OVR, quiet=QUIET)
        satbandno = satbandno + 1


    # Process SWIR
    # satellite band number (HYPERION SWIR [482-628])
    satbandno = 482  # Band 77 of HYPERION is first for atmospheric correction
    for idx in range(77,224+1):
        # Geometrical conditions (HYPERION VNIR = 32, SWIR = 33)
        geom = 33
        b = "rad." + str(idx)
        b_out = b.replace("rad.", "ref.")
        # Extract min max radiance stats
        rinfo = Module('r.info', map="rad."+str(idx), flags="s", stdout_=PIPE)
        rinfo_out = parse_key_val(rinfo.outputs.stdout)
        # Extract mean values from Dict
        minimum = float(rinfo_out["min"])
        maximum = float(rinfo_out["max"])
        if (minimum < 0.0):
            minimum = 0.0
        param = []
        # geometrical conditions=HYPERION
        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, a satellite)
        param.append(str("%.2f" % round(sens_height,2))+"\n")
        # i th band of HYPERION (atcorr internal no)
        param.append(str(satbandno))
        prm = "param_HYPERION.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(input=b_out, output=HYPERIONDir+"_"+b_out+".tif", overwrite=OVR, quiet=QUIET)
        satbandno = satbandno + 1

Signatures preparation and launching i.spec.sam

# Define the Hyperion files available for processing
listf="ref.8 ref.9 ref.10 ref.11 ref.12 ref.13 ref.14 ref.15 ref.16 ref.17 ref.18 ref.19 ref.20 ref.21 ref.22 ref.23 ref.24 ref.25 ref.26 ref.27 ref.28 ref.29 ref.30 ref.31 ref.32 ref.33 ref.34 ref.35 ref.36 ref.37 ref.38 ref.39 ref.40 ref.41 ref.42 ref.43 ref.44 ref.45 ref.46 ref.47 ref.48 ref.49 ref.50 ref.51 ref.52 ref.53 ref.54 ref.55 ref.56 ref.57 ref.77 ref.78 ref.79 ref.80 ref.81 ref.82 ref.83 ref.84 ref.85 ref.86 ref.87 ref.88 ref.89 ref.90 ref.91 ref.92 ref.93 ref.94 ref.95 ref.96 ref.97 ref.98 ref.99 ref.100 ref.101 ref.102 ref.103 ref.104 ref.105 ref.106 ref.107 ref.108 ref.109 ref.110 ref.111 ref.112 ref.113 ref.114 ref.115 ref.116 ref.117 ref.118 ref.119 ref.120 ref.121 ref.122 ref.123 ref.124 ref.125 ref.126 ref.127 ref.128 ref.129 ref.130 ref.131 ref.132 ref.133 ref.134 ref.135 ref.136 ref.137 ref.138 ref.139 ref.140 ref.141 ref.142 ref.143 ref.144 ref.145 ref.146 ref.147 ref.148 ref.149 ref.150 ref.151 ref.152 ref.153 ref.154 ref.155 ref.156 ref.157"

# Clean up .sign files
rm -f *.sign
# Name the signature file for i.spec.sam
sigfname="KL.sig"
rm -f $sigfname

echo "# Signatures for Hyperion image of KL" >> $sigfname

# Water in river
classname="water1"
coords="215455.8579625377,383556.12608450593"
fname="$classname.sign"
# Clean up temp file
rm -f temp
for file in $listf
do
    r.what map=$file coordinates=$coords separator='comma' >> temp
done
# Define cols number in sig file
lno=$(cat temp | wc -l)
sed -i 's/\(.*\),,\(.*\)/\2/' temp
tr '\n' ' ' < temp > $fname
echo "# $classname" >> $sigfname

# Forest area
classname="forest1"
coords="199906.1669829222,331418.4535104364"
fname="$classname.sign"
# Clean up temp file
rm -f temp
for file in $listf
do
    r.what map=$file coordinates=$coords separator='comma' >> temp
done
sed -i 's/\(.*\),,\(.*\)/\2/' temp
tr '\n' ' ' < temp > $fname
echo "# $classname" >> $sigfname

# Define rows number in sig file
signo=$(ls *.sign | wc -l)

# Complete the sig file
echo "# " >> $sigfname
echo "Matrix $signo by $lno" >> $sigfname 
n=0
for file in *.sign
do
    echo "row$n: $(cat $file)" >> $sigfname
    n=$(echo "$n + 1" | bc)
done

# Define the imagery working group
# If already existing and loaded, then empty it
i.group -l group=temp > temp
val=$(cat temp | wc -l)
if [ $val -ne "" ]
then i.group -r group=temp file=temp
fi
# Load it with the new set from $listf
i.group group=temp input=$(echo $listf | sed '/\ /,/') 

# Run a first set of Spectral Angle convergence
i.spec.sam group=temp input=$sigfname result=s
Retrieved from "https://grasswiki.osgeo.org/w/index.php?title=Hyperion&oldid=27331"