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 approach.

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.

fmask.landsatangles.bilinearInterp(xMin, xMax, yMin, yMax, cornerVals, x, y)[source]

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.

fmask.landsatangles.findCorners(info, inputs, outputs, otherargs)[source]

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.

fmask.landsatangles.findImgCorners(img, imgInfo)[source]

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)


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.

Return the lat/long of the centre of the image as

(ctrLat, ctrLong)


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….

fmask.landsatangles.makeAngles(info, inputs, outputs, otherargs)[source]

Called from RIOS

Make 4-layer sun and satellite angles for the image block

fmask.landsatangles.makeAnglesImage(templateimg, outfile, nadirLine, extentSunAngles, satAzimuth, imgInfo)[source]

Make a single output image file of the sun and satellite angles for every pixel in the template image.


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]

fmask.landsatangles.sunAnglesForExtent(imgInfo, mtlInfo)[source]

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.

fmask.landsatangles.sunAnglesForPoints(latDeg, longDeg, hourGMT, jdp)[source]

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.