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 - untar *.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 = rinfo_out["min"]
        maximum = 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
Retrieved from "https://grasswiki.osgeo.org/w/index.php?title=Hyperion&oldid=27326"