Source code for fmask.landsatTOA

#!/usr/bin/env python
"""
Module that handles convertion of scaled radiance (DN) 
values from USGS to Top of Atmosphere (TOA) reflectance (\*1000).
"""
# This file is part of 'python-fmask' - a cloud masking module
# Copyright (C) 2015  Neil Flood
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 3
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
from __future__ import print_function, division

import numpy
from osgeo import gdal
from rios import applier, cuiprogress, fileinfo
from . import config

gdal.UseExceptions()

# Derived by Pete Bunting from 6S
LANDSAT8_ESUN = [1876.61, 1970.03, 1848.9, 1571.3, 967.66, 245.73, 82.03, 361.72]
# From Chander, G., Markham, B.L., Helder, D.L. (2008)
# Summary of current radiometric calibration coefficients for Landsat MSS, TM, ETM+, and EO-1 ALI sensors
# http://landsathandbook.gsfc.nasa.gov/pdfs/Landsat_Calibration_Summary_RSE.pdf
LANDSAT4_ESUN = [1983.0, 1795.0, 1539.0, 1028.0, 219.8, 83.49]
LANDSAT5_ESUN = [1983.0, 1796.0, 1536.0, 1031.0, 220.0, 83.44]
LANDSAT7_ESUN = [1997.0, 1812.0, 1533.0, 1039.0, 230.8, 84.90]
# Derived by RSC, using Landsat-9 RSR & Kurucz solar spectrum
LANDSAT9_ESUN = [1881.97, 1972.61, 1851.75, 1572.03, 969.37, 246.18, 82.11, 360.48]

ESUN_LOOKUP = {'LANDSAT_4': LANDSAT4_ESUN,
            'LANDSAT_5': LANDSAT5_ESUN,
            'LANDSAT_7': LANDSAT7_ESUN,
            'LANDSAT_8': LANDSAT8_ESUN,
            'LANDSAT_9': LANDSAT9_ESUN}

RADIANCE_MULT = 'RADIANCE_MULT_BAND_%d'
RADIANCE_ADD = 'RADIANCE_ADD_BAND_%d'

LMAX_KEY = 'LMAX_BAND%d'
LMIN_KEY = 'LMIN_BAND%d'
QCALMAX_KEY = 'QCALMAX_BAND%d'
QCALMIN_KEY = 'QCALMIN_BAND%d'


# band numbers in mtl file for gain and offset for reflective
BAND_NUM_DICT = {'LANDSAT_4': (1, 2, 3, 4, 5, 7), 
    'LANDSAT_5': (1, 2, 3, 4, 5, 7),
    'LANDSAT_7': (1, 2, 3, 4, 5, 7),
    'LANDSAT_8': (1, 2, 3, 4, 5, 6, 7, 9),
    'LANDSAT_9': (1, 2, 3, 4, 5, 6, 7, 9)}


[docs]def readGainsOffsets(mtlInfo): """ Read the gains and offsets out of the .MTL file """ spaceCraft = mtlInfo['SPACECRAFT_ID'] nbands = len(BAND_NUM_DICT[spaceCraft]) gains = numpy.zeros(nbands) offsets = numpy.zeros(nbands) if (RADIANCE_MULT % 1) in mtlInfo: # Newer format for MTL file for idx, band in enumerate(BAND_NUM_DICT[spaceCraft]): s = RADIANCE_MULT % band gain = float(mtlInfo[s]) gains[idx] = gain s = RADIANCE_ADD % band offset = float(mtlInfo[s]) offsets[idx] = offset else: # Old format, calculate gain and offset from min/max values for (idx, band) in enumerate(BAND_NUM_DICT[spaceCraft]): lMax = float(mtlInfo[LMAX_KEY % band]) lMin = float(mtlInfo[LMIN_KEY % band]) qcalMax = float(mtlInfo[QCALMAX_KEY % band]) qcalMin = float(mtlInfo[QCALMIN_KEY % band]) gain = (lMax - lMin) / (qcalMax - qcalMin) offset = lMin - qcalMin * gain gains[idx] = gain offsets[idx] = offset return gains, offsets
[docs]def earthSunDistance(date): """ Given a date in YYYYMMDD will compute the earth sun distance in astronomical units """ import datetime year = int(date[:4]) month = int(date[4:6]) day = int(date[6:]) d1 = datetime.datetime(year, month, day) d2 = datetime.datetime(year, 1, 1) # first day of year deltaT = d1 - d2 jday = deltaT.days + 1 # julian day of year. ds = (1.0 - 0.01673 * numpy.cos(0.9856 * (jday - 4) * numpy.pi / 180.0)) return ds
[docs]def riosTOA(info, inputs, outputs, otherinputs): """ Called from RIOS """ nbands = inputs.infile.shape[0] infile = inputs.infile.astype(numpy.float32) inIgnore = otherinputs.inNull if inIgnore is None: inIgnore = 0 cosSunZen = numpy.cos(inputs.angles[3] * otherinputs.anglesToRadians) nullMask = (inputs.infile == inIgnore).any(axis=0) toaRefList = [] for band in range(nbands): rtoa = infile[band] * otherinputs.gains[band] + otherinputs.offsets[band] p = numpy.pi * rtoa * otherinputs.earthSunDistanceSq / (otherinputs.esun[band] * cosSunZen) # clip to a sensible range numpy.clip(p, 0.0, 2.0, out=p) toaRefList.append(p) out = numpy.array(toaRefList) * 10000.0 # convert to int16 outputs.outfile = out.astype(numpy.int16) # Mask out where input is null for i in range(len(outputs.outfile)): outputs.outfile[i][nullMask] = otherinputs.outNull
[docs]def makeTOAReflectance(infile, mtlFile, anglesfile, outfile): """ Main routine - does the calculation The eqn for TOA reflectance, p, is p = pi * L * d^2 / E * cos(theta) d = earthSunDistance(date) L = image pixel (radiance) E = exoatmospheric irradiance for the band, and theta = solar zenith angle. Assumes infile is radiance values in DN from USGS. mtlFile is the .mtl file. outfile will be created in the default format that RIOS is configured to use and will be top of atmosphere reflectance values *10000. Also assumes that the angles image file is scaled as radians*100, and has layers for satAzimuth, satZenith, sunAzimuth, sunZenith, in that order. """ mtlInfo = config.readMTLFile(mtlFile) spaceCraft = mtlInfo['SPACECRAFT_ID'] date = mtlInfo['DATE_ACQUIRED'] date = date.replace('-', '') inputs = applier.FilenameAssociations() inputs.infile = infile inputs.angles = anglesfile outputs = applier.FilenameAssociations() outputs.outfile = outfile otherinputs = applier.OtherInputs() otherinputs.earthSunDistance = earthSunDistance(date) otherinputs.earthSunDistanceSq = otherinputs.earthSunDistance * otherinputs.earthSunDistance otherinputs.esun = ESUN_LOOKUP[spaceCraft] gains, offsets = readGainsOffsets(mtlInfo) otherinputs.gains = gains otherinputs.offsets = offsets otherinputs.anglesToRadians = 0.01 otherinputs.outNull = 32767 imginfo = fileinfo.ImageInfo(infile) otherinputs.inNull = imginfo.nodataval[0] controls = applier.ApplierControls() controls.progress = cuiprogress.GDALProgressBar() controls.setStatsIgnore(otherinputs.outNull) controls.setCalcStats(False) applier.apply(riosTOA, inputs, outputs, otherinputs, controls=controls) # Explicitly set the null value in the output ds = gdal.Open(outfile, gdal.GA_Update) for i in range(ds.RasterCount): ds.GetRasterBand(i + 1).SetNoDataValue(otherinputs.outNull)