Hyperion: Difference between revisions
Jump to navigation
Jump to search
| Line 368: | Line 368: | ||
<source lang="shell"> | <source lang="shell"> | ||
# Define the Hyperion files available for processing | # 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" | 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" | ||
| Line 393: | Line 392: | ||
lno=$(cat temp | wc -l) | lno=$(cat temp | wc -l) | ||
sed -i 's/\(.*\),,\(.*\)/\2/' temp | sed -i 's/\(.*\),,\(.*\)/\2/' temp | ||
tr '\n' ' ' < | rm -f temp1 | ||
for number in $(cat temp) | |||
do | |||
printf "%.8f " $number >> temp1 | |||
done | |||
tr '\n' ' ' < temp1 > $fname | |||
echo "# $classname" >> $sigfname | echo "# $classname" >> $sigfname | ||
| Line 407: | Line 411: | ||
done | done | ||
sed -i 's/\(.*\),,\(.*\)/\2/' temp | sed -i 's/\(.*\),,\(.*\)/\2/' temp | ||
tr '\n' ' ' < | rm -f temp1 | ||
for number in $(cat temp) | |||
do | |||
printf "%.8f " $number >> temp1 | |||
done | |||
tr '\n' ' ' < temp1 > $fname | |||
echo "# $classname" >> $sigfname | echo "# $classname" >> $sigfname | ||
| Line 415: | Line 424: | ||
# Complete the sig file | # Complete the sig file | ||
echo "# " >> $sigfname | echo "# " >> $sigfname | ||
echo "Matrix $signo by $lno" >> $sigfname | echo "Matrix: $signo by $lno" >> $sigfname | ||
n=0 | n=0 | ||
for file in *.sign | for file in *.sign | ||
| Line 425: | Line 434: | ||
# Define the imagery working group | # Define the imagery working group | ||
listmaps=$(echo $listf | sed 's/\ /,/g') | listmaps=$(echo $listf | sed 's/\ /,/g') | ||
i.group group=temp input=$listmaps | i.group group=temp input=$listmaps --quiet | ||
# Run a first set of Spectral Angle convergence | # Run a first set of Spectral Angle convergence | ||
i.spec.sam group=temp input=$sigfname result=s | i.spec.sam group=temp input=$sigfname result=s | ||
</source> | </source> | ||
Revision as of 11:50, 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.
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).
| 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
rm -f temp1
for number in $(cat temp)
do
printf "%.8f " $number >> temp1
done
tr '\n' ' ' < temp1 > $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
rm -f temp1
for number in $(cat temp)
do
printf "%.8f " $number >> temp1
done
tr '\n' ' ' < temp1 > $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
listmaps=$(echo $listf | sed 's/\ /,/g')
i.group group=temp input=$listmaps --quiet
# Run a first set of Spectral Angle convergence
i.spec.sam group=temp input=$sigfname result=s