Hyperion: Difference between revisions
m (→Data Download) |
|||
Line 53: | Line 53: | ||
Of which ''EO1H1260582015073110T1_MTL_L1GST.TXT'' is a Landsat type metadata header file (''MTL type'') well known to Landsat users. | Of which ''EO1H1260582015073110T1_MTL_L1GST.TXT'' is a Landsat type metadata header file (''MTL type'') well known to Landsat users. | ||
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. | === 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 == | == Import in GRASS, DN2Rad & Rad2Ref == |
Revision as of 06:57, 2 December 2023
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.
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).
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 = 0
# 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 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 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')
exit()
# 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 = 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
# 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