Source code for fmask.landsatangles

#!/usr/bin/env python
# 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
# 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.
Functions relating to estimating the per-pixel sun and satellite angles for 
a given Landsat image. These are rough estimates, using the generic
characteristics of the Landsat 5 platform, and are not particularly accurate,
but good enough for the current purposes. 

Historically, the USGS have not supplied satellite zenith/azimuth angles, and have only 
supplied scene-centre values for sun zenith/azimuth angles. Since the satellite
view geometry is important in correctly tracking a shadow when matching shadows
to their respective clouds, the Fmask algorithm requires good estimates of all these
angles. The routines contained here are used to derive per-pixel estimates of 
these angles. 

As of mid-2016, the USGS are planning to supply sufficient information to calculate
these angles directly from orbit ephemeris data. When that comes about, it seems likely
that the need for the routines here will diminish, but any data downloaded from USGS
prior to then will still require this approach, as the associated angle metadata will 
not be present. 

The core Fmask code in this package is adaptable enough to be configured for either 

The general approach for satellite angles is to estimate the nadir line by running it
down the middle of the image data area. The satellite azimuth is assumed to be
at right angles to this nadir line, which is only roughly correct. For the whisk-broom 
sensors on Landsat-5 and Landsat-7, this angle is not 90 degrees, but is affected by 
earth rotation and is latitude dependent. For Landsat-8, the scan line is at 
right angles, due to the compensation for earth rotation, but the push-broom is 
made up of sub-modules which point in slightly different directions, giving 
slightly different satellite azimuths along the scan line. None of these effects
are included in the current estimates. The satellite zenith is estimated based on the
nadir point, the scan-line, and the assumed satellite altitude, and includes the
appropriate allowance for earth curvature. 

Because this works by searching the imagery for the non-null area, and assumes that 
this represents a full-swath image, it would not work for a subset of a full image. 

The sun angles are approximated using the algorithm found in the Fortran code with
6S (Second Simulation of the Satellite Signal in the Solar Spectrum). The subroutine
in question is the POSSOL() routine. I translated the Fortran code into Python for
inclusion here. 

from __future__ import print_function, division

import datetime

import numpy
from osgeo import osr

from rios import applier
from rios import fileinfo

[docs]def findImgCorners(img, imgInfo): """ Find the corners of the data within the given template image Return a numpy array of (x, y) coordinates. The array has 2 columns, for X and Y. Each row is a corner, in the order: top-left, top-right, bottom-left, bottom-right. Uses RIOS to pass through the image searching for non-null data, and find the extremes. Assumes we are working with a full-swathe Landsat image. Each list element is a numpy array of (x, y) """ infiles = applier.FilenameAssociations() outfiles = applier.FilenameAssociations() otherargs = applier.OtherInputs() infiles.img = img = None = None = None = None otherargs.nullVal = imgInfo.nodataval[0] if otherargs.nullVal is None: otherargs.nullVal = 0 applier.apply(findCorners, infiles, outfiles, otherargs) corners = numpy.array([,,,, ]) return corners
[docs]def findCorners(info, inputs, outputs, otherargs): """ Called from RIOS Checks non-null area of image block. Finds extremes, records coords of extremes against those already in otherargs. Note that the logic is very specific to the orientation of the usual Landsat descending pass imagery. The same logic should not be applied to swathes oriented in other, e.g. for other satellites. """ (xblock, yblock) = info.getBlockCoordArrays() nonnull = (inputs.img != otherargs.nullVal).all(axis=0) xNonnull = xblock[nonnull] yNonnull = yblock[nonnull] if len(xNonnull) > 0: topNdx = numpy.argmax(yNonnull) topXY = (xNonnull[topNdx], yNonnull[topNdx]) leftNdx = numpy.argmin(xNonnull) leftXY = (xNonnull[leftNdx], yNonnull[leftNdx]) bottomNdx = numpy.argmin(yNonnull) bottomXY = (xNonnull[bottomNdx], yNonnull[bottomNdx]) rightNdx = numpy.argmax(xNonnull) rightXY = (xNonnull[rightNdx], yNonnull[rightNdx]) # If these are more extreme than those already in otherargs, replace them if is None or topXY[1] >[1]: = topXY if is None or rightXY[0] >[0]: = rightXY if is None or leftXY[0] <[0]: = leftXY if is None or bottomXY[1] <[1]: = bottomXY
[docs]def findNadirLine(corners): """ Return the equation of the nadir line, from the given corners of the swathe. Returns a numpy array of [b, m], for the equation y = mx + b giving the y coordinate of the nadir as a function of the x coordinate. """ # Find the top and bottom mid-points. topMid = (corners[0] + corners[1]) / 2.0 bottomMid = (corners[2] + corners[3]) / 2.0 slope = (topMid[1] - bottomMid[1]) / (topMid[0] - bottomMid[0]) intercept = topMid[1] - slope * topMid[0] coeffs = numpy.array([intercept, slope]) return coeffs
[docs]def satAzLeftRight(nadirLine): """ Calculate the satellite azimuth for the left and right sides of the nadir line. Assume that the satellite azimuth vector is at right angles to the nadir line (which is not really true, but reasonably close), and that there are only two possibilities, as a pixel is either to the left or to the right of the nadir line. Return a numpy array of [satAzLeft, satAzRight], with angles in radians, in the range [-pi, pi] """ slope = nadirLine[1] # Slope of a line perpendicular to the nadir slopePerp = -1 / slope # Azimuth for pixels to the left of the line azimuthLeft = numpy.pi / 2.0 - numpy.arctan(slopePerp) # Azimuth for pixels to the right is directly opposite azimuthRight = azimuthLeft - numpy.pi return numpy.array([azimuthLeft, azimuthRight])
[docs]def localRadius(latitude): """ Calculate a local radius of curvature, for the given geodetic latitude. This approximates the earth curvature at this latitude. The given latitude is in degrees. This is probably overkill, given some of the other approximations I am making.... """ latRadians = numpy.radians(latitude) # Earth axis lengths a = osr.SRS_WGS84_SEMIMAJOR invFlat = osr.SRS_WGS84_INVFLATTENING f = 1 / invFlat eSqd = 2*f - f**2 # Radius of curvature R = a / numpy.sqrt(1 - eSqd * numpy.sin(latRadians)**2) return R
[docs]def sunAnglesForExtent(imgInfo, mtlInfo): """ Return array of sun azimuth and zenith for each of the corners of the image extent. Note that this is the raster extent, not the corners of the swathe. The algorithm used here has been copied from the 6S possol() subroutine. The Fortran code I copied it from was .... up to the usual standard in 6S. So, the notation is not always clear. """ cornerLatLong = imgInfo.getCorners(outEPSG=4326) (ul_long, ul_lat, ur_long, ur_lat, lr_long, lr_lat, ll_long, ll_lat) = cornerLatLong pts = numpy.array([ [ul_long, ul_lat], [ur_long, ur_lat], [ll_long, ll_lat], [lr_long, lr_lat] ]) longDeg = pts[:, 0] latDeg = pts[:, 1] # Date/time in UTC dateStr = mtlInfo['DATE_ACQUIRED'] timeStr = mtlInfo['SCENE_CENTER_TIME'].replace('Z', '') ymd = [int(i) for i in dateStr.split('-')] dateObj =[0], ymd[1], ymd[2]) julianDay = (dateObj -[0], 1, 1)).days + 1 juldayYearEnd = ([0], 12, 31) -[0], 1, 1)).days + 1 # Julian day as a proportion of the year jdp = julianDay / juldayYearEnd # Hour in UTC hms = [float(x) for x in timeStr.split(':')] hourGMT = hms[0] + hms[1] / 60.0 + hms[2] / 3600.0 (sunAz, sunZen) = sunAnglesForPoints(latDeg, longDeg, hourGMT, jdp) sunAngles = numpy.vstack((sunAz, sunZen)).T return sunAngles
[docs]def sunAnglesForPoints(latDeg, longDeg, hourGMT, jdp): """ Calculate sun azimuth and zenith for the given location(s), for the given time of year. jdp is the julian day as a proportion, ranging from 0 to 1, where Jan 1 is 1.0/365 and Dec 31 is 1.0. Location is given in latitude/longitude, in degrees, and can be arrays to calculate for multiple locations. hourGMT is a decimal hour number giving the time of day in GMT (i.e. UTC). Return a tuple of (sunAz, sunZen). If latDeg and longDeg are arrays, then returned values will be arrays of the same shape. """ latRad = numpy.radians(latDeg) # Express jdp in radians jdpr = jdp * 2 * numpy.pi # Now work out the solar position. This is copied from the 6S code, but # is also documented in the 6S manual. The notation a = numpy.array([0.000075, 0.001868, 0.032077, 0.014615, 0.040849]) meanSolarTime = hourGMT + longDeg / 15.0 localSolarDiff = (a[0] + a[1] * numpy.cos(jdpr) - a[2] * numpy.sin(jdpr) - a[3] * numpy.cos(2*jdpr) - a[4] * numpy.sin(2*jdpr)) * 12 * 60 / numpy.pi trueSolarTime = meanSolarTime + localSolarDiff / 60 - 12.0 # Hour as an angle ah = trueSolarTime * numpy.radians(15) b = numpy.array([0.006918, 0.399912, 0.070257, 0.006758, 0.000907, 0.002697, 0.001480]) delta = (b[0] - b[1]*numpy.cos(jdpr) + b[2]*numpy.sin(jdpr) - b[3]*numpy.cos(2.*jdpr) + b[4]*numpy.sin(2.*jdpr) - b[5]*numpy.cos(3.*jdpr) + b[6]*numpy.sin(3.*jdpr)) cosSunZen = (numpy.sin(latRad) * numpy.sin(delta) + numpy.cos(latRad) * numpy.cos(delta) * numpy.cos(ah)) sunZen = numpy.arccos(cosSunZen) # sun azimuth from south, turning west (yeah, I know, weird, isn't it....) sinSunAzSW = numpy.cos(delta) * numpy.sin(ah) / numpy.sin(sunZen) sinSunAzSW = sinSunAzSW.clip(-1.0, 1.0) # This next bit seems to be to get the azimuth in the correct quadrant. I do # not really understand it. cosSunAzSW = (-numpy.cos(latRad) * numpy.sin(delta) + numpy.sin(latRad) * numpy.cos(delta) * numpy.cos(ah)) / numpy.sin(sunZen) sunAzSW = numpy.arcsin(sinSunAzSW) sunAzSW = numpy.where(cosSunAzSW <= 0, numpy.pi - sunAzSW, sunAzSW) sunAzSW = numpy.where((cosSunAzSW > 0) & (sinSunAzSW <= 0), 2 * numpy.pi + sunAzSW, sunAzSW) # Now convert to azimuth from north, turning east, as is usual convention sunAz = sunAzSW + numpy.pi # Keep within [0, 2pi] range sunAz = numpy.where(sunAz > 2 * numpy.pi, sunAz - 2 * numpy.pi, sunAz) return (sunAz, sunZen)
[docs]def makeAnglesImage(templateimg, outfile, nadirLine, extentSunAngles, satAzimuth, imgInfo): """ Make a single output image file of the sun and satellite angles for every pixel in the template image. """ imgInfo = fileinfo.ImageInfo(templateimg) infiles = applier.FilenameAssociations() outfiles = applier.FilenameAssociations() otherargs = applier.OtherInputs() controls = applier.ApplierControls() infiles.img = templateimg outfiles.angles = outfile (ctrLat, ctrLong) = getCtrLatLong(imgInfo) otherargs.R = localRadius(ctrLat) otherargs.nadirLine = nadirLine otherargs.xMin = imgInfo.xMin otherargs.xMax = imgInfo.xMax otherargs.yMin = imgInfo.yMin otherargs.yMax = imgInfo.yMax otherargs.extentSunAngles = extentSunAngles otherargs.satAltitude = 705000 # Landsat nominal altitude in metres otherargs.satAzimuth = satAzimuth otherargs.radianScale = 100 # Store pixel values as (radians * radianScale) controls.setStatsIgnore(500) applier.apply(makeAngles, infiles, outfiles, otherargs, controls=controls)
[docs]def makeAngles(info, inputs, outputs, otherargs): """ Called from RIOS Make 4-layer sun and satellite angles for the image block """ (xblock, yblock) = info.getBlockCoordArrays() # Nadir line coefficients of y=mx+b (b, m) = otherargs.nadirLine # Distance of each pixel from the nadir line dist = numpy.absolute((m * xblock - yblock + b) / numpy.sqrt(m**2 + 1)) # Zenith angle assuming a flat earth satZenith = numpy.arctan(dist / otherargs.satAltitude) # Adjust satZenith for earth curvature. This is a very simple approximation, but # the adjustment is less than one degree anyway, so this is accurate enough. curveAngle = numpy.arctan(dist / otherargs.R) satZenith += curveAngle # Work out whether we are left or right of the nadir line isLeft = (yblock - (m * xblock + b)) > 0 (satAzimuthLeft, satAzimuthRight) = otherargs.satAzimuth satAzimuth = numpy.where(isLeft, satAzimuthLeft, satAzimuthRight) # Interpolate the sun angles from those calculated at the corners of the whole raster extent (xMin, xMax, yMin, yMax) = (otherargs.xMin, otherargs.xMax, otherargs.yMin, otherargs.yMax) sunAzimuth = bilinearInterp(xMin, xMax, yMin, yMax, otherargs.extentSunAngles[:, 0], xblock, yblock) sunZenith = bilinearInterp(xMin, xMax, yMin, yMax, otherargs.extentSunAngles[:, 1], xblock, yblock) angleStack = numpy.array([satAzimuth, satZenith, sunAzimuth, sunZenith]) angleStackDN = angleStack * otherargs.radianScale outputs.angles = numpy.round(angleStackDN).astype(numpy.int16)
[docs]def bilinearInterp(xMin, xMax, yMin, yMax, cornerVals, x, y): """ Evaluate the given value on a grid of (x, y) points. The exact value is given on a set of corner points (top-left, top-right, bottom-left, bottom-right). The corner locations are implied from xMin, xMax, yMin, yMax. """ p = (y - yMin) / (yMax - yMin) q = (x - xMin) / (xMax - xMin) # Give the known corner values some simple names (tl, tr, bl, br) = cornerVals # Calculate the interpolated values vals = tr * p * q + tl * p * (1-q) + br * (1-p) * q + bl * (1-p) * (1-q) return vals
[docs]def getCtrLatLong(imgInfo): """ Return the lat/long of the centre of the image as (ctrLat, ctrLong) """ cornerLatLong = imgInfo.getCorners(outEPSG=4326) (ul_long, ul_lat, ur_long, ur_lat, lr_long, lr_lat, ll_long, ll_lat) = cornerLatLong ctrLat = numpy.array([ul_lat, ur_lat, lr_lat, ll_lat]).mean() ctrLong = numpy.array([ul_long, ur_long, lr_long, ll_long]).mean() return (ctrLat, ctrLong)