Remote Sensing Applications: Society and Environment 6 (2017) 1–14

Contents lists available at ScienceDirect

Remote Sensing Applications: Society and Environment

journal homepage: www.elsevier.com/locate/rsase

A standardised Landsat time series (1973–2016) of forest leaf area index MARK
using pseudoinvariant features and spectral vegetation index isolines and a
catchment hydrology application
⁎
Craig Macfarlanea, , Andrew H. Griggb, Matthew I. Dawsb
a
CSIRO, 147 Underwood Avenue, Floreat, WA 6014, Australia
b
Environmental Department, Alcoa of Australia, PO Box 172, Pinjarra, WA 6208, Australia

A R T I C L E I N F O A BS T RAC T

Keywords: Remotely-sensed imagery from Landsat dating back to 1972 is now available from the United States Geological
Multispectral Scanner (MSS) Survey free-of-charge as a level 1 terrain corrected product. However, to take full advantage of the time series
Thematic Mapper (TM, ETM+) requires that the earlier multispectral scanner (MSS) imagery be integrated with the later thematic mapper
Operational Land Imager (OLI) (TM), enhanced thematic mapper (ETM+) and operational land imager (OLI) imagery. Here we describe a
Atmospheric correction
simple, generic approach to processing Landsat scenes to develop a standardised time series of leaf area index, a
Bare soil line
key biophysical attribute of forests, for a study region in the northern jarrah forest of south-western Australia.
Simple ratio
NDVI The five-step approach utilised the main features in near infra-red (NIR)-red space of a dark point, soil line and
Leaf area index a radiating set of ratio index isolines. Firstly, digital numbers were converted to top of atmosphere reflectance.
Cover photography Secondly, the red and NIR bands of all scenes were atmospheric corrected using long-established deep water
Hemispherical photography supply reservoirs within the study region as invariant ‘dark objects’. Thirdly, adjustment of the NIR band such
Streamflow that the soil line was consistently located on the 1:1 line in NIR-red space across scenes both from the same
Jarrah forest sensor and from different sensors, taking advantage of a > 40 year record of bauxite mining operations to
identify pseudo-invariant bare-ground targets. Fourthly, an inter-sensor calibration of spectral vegetation
indices (SVIs) to ensure that SVI isolines from different sensors were consistent. Finally, application of a
relationship between SVIs and an extensive multi-year dataset of ground-based LAI measurements to generate a
time-series of leaf area index (LAI). We demonstrate a simple application of the LAI time-series by investigating
the relationship between rainfall, forest LAI as affected by different land use activities, and streamflow in the
study region in the south west of Australia. The time-series could be used to improve prediction and modeling of
hydrological changes in the region resulting from changes in catchment forest cover.

1. Introduction these are impractical at regional scales. The often strong correlations
between spectral vegetation indices (SVIs) and forest dynamics (Cohen
Remote sensing of leaf area index (LAI) in forest ecosystems is and Goward, 2004) make it possible to assess the potential impacts of
important because the forest canopy is the main interface through large scale deforestation and reforestation on the Earth's carbon cycle
which carbon is sequestered into biomass via photosynthesis, through (Canadell and Raupach, 2008), the likely impacts that on-going climate
which water is lost from forests to the atmosphere via evaporation, and change is in-turn predicted to have on the world's forests (Aitken et al.,
through which vegetation alters the energy exchange between the land 2008), and other large scale anthropogenic impacts on forests.
surface and the atmosphere. LAI is a key determinant of the produc- Remotely-sensed imagery from the various Landsat missions dates
tivity and water use of vegetation (Watson, 1958; Lieth, 1975; Waring, back to 1972. In addition to its long historical archive, the Landsat time
1983; Gholz et al., 1990; Baldocchi and Harley, 1995; Sellers et al., series (LTS) is well suited to time series analysis owing to the relatively
1997; Ryu et al., 2011) and is also a key input to process-based fine spatial resolution of its multispectral channels (circa 30 m from
ecosystem models (e.g., Running and Gower, 1991; Landsberg and Landsat 4 onwards; Markham and Helder, 2012). Furthermore, all new
Waring, 1997; Wang et al., 2011). Whilst field methods exist to and archived Landsat imagery is now available from the United States
measure LAI at local scales (Breda, 2003; Macfarlane et al., 2007b), Geological Survey (USGS) free-of-charge as a level 1 terrain corrected

⁎
Corresponding author.
E-mail address: craig.macfarlane@csiro.au (C. Macfarlane).

http://dx.doi.org/10.1016/j.rsase.2017.01.006
Received 30 August 2016; Accepted 22 January 2017
Available online 23 January 2017
2352-9385/ Crown Copyright © 2017 Published by Elsevier B.V. All rights reserved.
C. Macfarlane et al. Remote Sensing Applications: Society and Environment 6 (2017) 1–14

product (Woodcock et al., 2008; Wulder et al., 2012). However, to take studies (Powell et al., 2008; Gómez et al., 2011; Pflugmacher et al.,
full advantage of the LTS requires that the earlier multispectral scanner 2012), we sought to incorporate MSS imagery with later Landsat
(MSS) imagery be integrated with the later thematic mapper (TM), imagery using an approach based on atmospheric correction of the top-
enhanced thematic mapper (ETM+) and, since 2014, operational land of-atmosphere reflectance, adjustment of the bare soil line and inter-
imager (OLI) imagery. The MSS imagery supplied by the USGS since sensor calibration of vegetation index isolines. The time-series of LAI
2010 has improved radiometric calibrations and improved geometric was used to investigate the relationship between rainfall, forest LAI as
corrections (Markham and Helder, 2012) thus improving its consis- affected by land use activities, and streamflow. The northern jarrah
tency with later Landsat sensors; however, the earlier MSS sensor has forest is an extensive tract of relatively intact forest cover that provides
lesser spatial resolution (circa 60 m) than the later sensors (Markham a range of land uses including conservation, timber production, bauxite
and Helder, 2012) and has markedly different bandwidth and band mining and water supply. Catchments within the forest are an
positioning to the later sensors, especially in the near infra-red (NIR). important source of fresh drinking water, and there has been a long
Markham and Helder (2012, p.18) noted that “band 4 is much broader interest in the effect on catchment flows and water quality due to
with the MSS sensor, and band 3 MSS really has no counterpart with changes in forest cover or catchment LAI (Stoneman and Schofield,
respect to TM”. Given the importance of the NIR band to studies of 1989; Ruprecht and Stoneman, 1993). We further tested whether the
ecology and vegetation (Cohen and Goward, 2004), this is a significant response of streamflow to change in LAI in catchments subject to
obstacle to integration of the MSS imagery with the rest of the LTS, and bauxite mining can be distinguished from the response in catchments
few long time-series studies have utilised MSS imagery (Hostert et al., subject to other land uses (Bari and Ruprecht, 2003).
2003; Powell et al., 2008; Gómez et al., 2011; Pflugmacher et al., 2012).
Indeed, Gallo and Daughtry (1987, p.440) stated “Differences between 2. Materials and methods
sensor systems in visible and near-IR wavebands used to compute
vegetation indices may be too great, however, to permit direct 2.1. Study region
comparison” and Markham and Helder (2012, p.18) stated “it must
be emphasized that MSS and TM cannot be expected to give the same The study region, encompassing 27,000 km2 in the south-west of
radiometric response to targets with varying spectral signatures”. Australia, is centred on the northern jarrah forest with the city of Perth
Galvao et al. (1999) found that NDVI values were affected by the (31.9522°S, 115.8589°E, population > 2 million people), on its wes-
proximity of the red and NIR bands to the spectral interval of the red tern edge (Fig. 1). The distinct western forest boundary aligns with a
edge (690–750 nm), which is likely to be a particular problem for 300 m elevational shift associated with the Darling Plateau escarpment,
sensors with broad spectra such as MSS. Previous studies have also responsible for a strong rainfall gradient which declines from west to
found that relative calibrations of different sensors may be specific to east (Schofield et al., 1989). The higher rainfall ( > 1100 mm long-term
different types of land cover, although this is more the case for annual average) parts of the forest on this western boundary have
individual bands than for derived indices such as NDVI (Yoshioka, historically provided the majority of Perth's water supply, and several
2003; Miura et al., 2006). dams and reservoirs have been constructed since the early 1900's
The tasselled cap transformation (TCT) has been a popular choice (Fig. 1). South-western Australia has experienced a significant decline
for previous studies that utilised MSS imagery in LTS (Powell et al., in rainfall since the mid-1970's (Bates et al., 2008) and many once
2008; Gómez et al., 2011; Pflugmacher et al., 2012) for the purposes of perennial streams in the northern jarrah forest have now become
studying forest dynamics. The TCT utilises 4–6 bands to calculate three seasonal (Petrone et al., 2010), with further declines in both rainfall
components, two of which, the greenness and brightness components, and streamflow expected (Silberstein et al., 2012). The forest is
form a ‘vegetation plane’ (Kauth et al., 1976). Indices such as the dominated by the overstorey trees jarrah and marri (Corymbia
‘tasselled cap angle’ can condense this information into a single value calophylla) with a diverse sclerophyllous understorey of small to
that can be related to forest structure (e.g. Gómez et al., 2011). medium shrubs and ground layer species (Bell and Heddle, 1989).
However, the choice of spectral vegetation index (SVI) depends upon Soils throughout the study region are relatively uniform and deeply ( >
the application. SVIs such as normalized-difference vegetation index 30 m) weathered comprised of lateritic gravels, sands and loams
(NDVI; Rouse et al., 1974) and simple ratio (SR, Jordan, 1969), which overlying yellow to white mottled and pallid zone clays and basement
are based on the much larger reflectance of NIR than red by vegetation, rock of predominantly Archean granite.
are more commonly used for prediction of LAI from remote sensing Mining for bauxite is a major land use within the jarrah forest,
(e.g., Baret and Guyot, 1991; Coops et al., 1997; Peddle et al., 1999; commencing in the early 1960's and expanding rapidly through the
Turner et al., 1999; Qi et al., 2000; Gascon et al., 2004; Soudani et al., 1970's (Havel, 1989). The soils contain extensive but discontinuous
2006; Boer et al., 2008). In NIR-red space there is a principal axis of and relatively shallow (3–5 m) deposits of bauxite (Hickman et al.,
soil spectral variation, or “soil line”, extending from the origin, or “dark 1992) which have been mined and progressively rehabilitated back to
point”, with increasing brightness (Baret et al., 1993; Huete, 1988) and forest cover since the 1960's (Gardner and Bell, 2007; Koch, 2007).
ratio SVIs such as the SR or NDVI are represented by isolines of While clearing and subsequent restoration during mining are similar to
increasing slope (NIR/red) extending from the origin above the soil other disturbances such as logging or thinning in terms of changes in
line. These are also isolines of constant vegetation amount, provided forest cover, excavation of the upper soil profile is not, and catchments
that background soil characteristics are also constant (Huete, 1988). with significant proportions of mining may therefore respond hydro-
Hence, in order to construct a standardised time series of LAI from SR logically in a different way to other land uses.
or NDVI it is necessary to standardize the location of the dark point
and soil line, and also ensure that isolines of NDVI and SR across 2.2. Approach, scene selection and processing
scenes represent the same variation of vegetation amount. Adjustment
of the dark point is known as atmospheric correction because it 2.2.1. Approach
corrects for atmospheric scattering or haze, which alters the reflectance We adopted a five-step approach to develop a standardised time-
sensed by the satellite (Chavez, 1988; Song et al., 2001; Chang et al., series of LAI from Landsat imagery using the main features in NIR-red
2008). space of a dark point, soil line and a radiating set of ratio index isolines
In this study, we aimed to develop a standardized annual time- as shown in Fig. 2:
series of LAI for the northern jarrah (Eucalyptus marginata) forest of
south-western Australia using the 44 year time series of level 1 terrain 1. We first converted the digital numbers to top of atmosphere
corrected Landsat imagery sourced from the USGS. Unlike previous reflectance based on the supplied calibration parameters.

2
C. Macfarlane et al. Remote Sensing Applications: Society and Environment 6 (2017) 1–14

Fig. 1. Study region in south-western Australia (inset). The dark green areas within the satellite imagery are forested while the paler areas are mainly agricultural and, around major
centres, urban or industrial land use. The black rectangle indicates the study region within the two mosaicked Landsat 4–8 scenes (light blue). The locations of four deep water supply
reservoirs used both as ‘dark points’ for atmospheric correction and for pixel registration are indicated in blue, and the locations of 30 gauged catchments (Table 2) used in the
development of the streamflow model are indicated in red. Detail of the reservoir and catchment locations is shown on the right. This map was created using the World Imagery basemap
within ArcGIS® software by Esri. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

2. We applied a simple dark object subtraction (DOS; Chavez, 1988) both from the same sensor and from different sensors. We took
atmospheric correction to the red and NIR bands of all scenes using advantage of the known locations of bare ground based on a > 40
long-established water supply reservoirs within the study region year record of bauxite mining operations in the study region (Alcoa,
(Fig. 1) that contain deep, clear water as invariant targets. unpublished data) to identify pseudo-invariant bare-ground targets.
3. We then adjusted the NIR band such that the soil line was 4. We then tested whether calculated NDVI and SR isolines across
consistently located on the 1:1 line in NIR-red space across scenes different sensors represented the same variation of vegetation

3
C. Macfarlane et al. Remote Sensing Applications: Society and Environment 6 (2017) 1–14

the series owing either to excessive ( > 10%) cloud cover over the study
region or unavailability for other reasons. The summer months are also
when incoming radiance is highest, allowing for greater reflectance-
based discrimination, which is important for detecting change
(McVicar et al., 2002). For Landsat 1–3, most of the study region
was obtained from path/row 120/82 and 119/83 a day apart, but in
cases where scenes with widely different dates (and in some cases from
different sensors) were used to avoid cloud, scene brightness was
adjusted based on a comparison of overlapping pixels. For Landsat 4 –
8, we selected scenes along path 112 and rows 82 and 83 that were
acquired on the same date. Individual scenes were mosaicked and then
cropped to the region of interest (Fig. 1). The cropped mosaics were
converted to TOA reflectance, corrected for solar elevation angle, as
described in the Landsat User Handbooks (US Geological Survey,
2012; US Geological Survey, 2015): Landsat 8 DNs were converted
directly to TOA reflectance while the DNs of other sensors were first
converted to radiance and then to reflectance. The mosaics were saved
as 16 bit TIFFs such that zero reflectance was zero and a reflectance of
one was represented as 65535.
Fig. 2. Schematic diagram of NIR-red space for the calculation of SR and NDVI. The soil Red and NIR bands used were bands 5 and 6 for Landsat 1–3 MSS
line is the set of points at which SR=1 and NDVI=0. sensors, bands 2 and 3 for Landsat 4–5 MSS sensors, bands 3 and 4 for
Landsat 4–5 TM and Landsat 7 ETM+ sensors and bands 4 and 5 for
amount, by comparing SVIs from different sensors at the same the Landsat 8 OLI sensor. We selected band 3 from MSS 4–5 to
location, at either the same time or eight days apart. We developed represent NIR (equivalent to band 6 from MSS 1–3) based on the more
inter-sensor calibrations of SVIs to ensure that SVI isolines from certain cross-calibration of this band compared to the longer wave-
different sensors were consistent. length band 4 (Markham and Helder, 2012). Hall et al. (1991) reported
5. Finally, we determined and applied a relationship between SVIs and that, post-launch, calibrations of Landsat 2–5 MSS only agreed to
ground-based LAI measurements to generate a time-series of LAI. within 12% of one another and tended to diverge further with age.
Recent cross-calibration (Markham and Helder, 2012) has reduced
2.2.2. Scene selection and registration uncertainty in the three visible MSS bands to ~4%, while the NIR band
During August 2014 to December 2015 Landsat scenes (level 1 still has much greater uncertainty (~12%).
terrain (L1T) corrected product) covering the region (Fig. 1) were Pflugmacher et al. (2012) found that, despite being terrain
obtained from the USGS EarthExplorer website (http://earthexplorer. corrected, the geometric accuracy of MSS was not as high as TM and
usgs.gov/) to provide an annual series from December 1972 to ETM+ and required further co-registration to a reference image. We
December 2015 (Table S1). The L1T product is radiometrically and geometrically co-registered all scenes in the series to a Landsat 8 scene
geometrically corrected based on inputs from the sensors and the from 10th January 2014, which we refer to as our reference scene.
spacecraft, as well as ground control points and digital elevation Band 8 (panchromatic with 15 m resolution) was used to evaluate the
models. The L1T image is presented in units of Digital Numbers absolute registration of the reference scene by aligning a point in the
(DNs) which can be rescaled to spectral radiance or top of atmosphere scene to a single precisely known ground location. Geometric correc-
(TOA) reflectance, based on calibration parameters supplied with the tions were automatically applied to other scenes as required using
imagery. MSS, TM and ETM+ scenes are supplied as 8 bit GeoTIFFs, intensity-based image registration functions (imregtform and imregis-
consistent with the 8 bit radiometric density of TM and ETM+ ter) in the MATLAB Image Processing Toolbox (MATLAB Release
(Table 1). The original radiometric density of MSS is only 7 bits but 2016a, The MathWorks, Inc., Natick, Massachusetts, United States),
MSS L1T scenes are rescaled to 8 bits. The radiometric density of OLI based on four long-established ( > 40 years) reservoirs from within the
is 12 bits and OLI products are supplied as 16 bit GeoTIFFs. study region (Fig. 1) which had high-contrast boundaries with sur-
Specifications of the various Landsat missions are given in Table 1. rounding vegetated areas. For each scene, the location of these
Where possible, scenes were selected within the months of reservoirs was compared to the reference scene to determine offsets
December to February (i.e., the austral summer) when the sun was for each scene, which were smaller than 30 m for all TM, ETM+ and
high overhead to minimise shadows and to minimise the variation in OLI scenes and smaller than 60 m for all MSS scenes.
solar angle. Scenes outside this period were necessary for four years in

Table 1
Spectral range (nm), dynamic range (W m2 sr−1 µm−1) and radiometric density (number of bits) for the red and NIR bands of Landsat missions 1–8. Specifications for Landsat 1 – 7
taken from Chander et al. (2009). Specifications for Landsat 8 taken from Morfitt et al. (2015).

Landsat mission Red band NIR band

nm W m2 sr−1 µm−1 # bits nm W m2 sr−1 µm−1 # bits

1 MSS 603–701 0–200 7-bit 694–800 0–176 7-bit

2 MSS 607–710 6–176 7-bit 697–802 6–152 7-bit
3 MSS 606–705 3–179 7-bit 693–793 3–149 7-bit
4 MSS 603–696 4–164 7-bit 701–813 5–142 7-bit
5 MSS 603–697 3–179 7-bit 704–814 5–148 7-bit
4 TM 624–693 −1.17–254 8-bit 776–905 −1.51–221 8-bit
5 TM 626–693 −1.17–264 8-bit 776–904 −1.51–221 8-bit
7 ETM+ 631–692 −5.0–234.4 8-bit 772–898 −5.1–241.1 8-bit
8 OLI 630–680 0.08–462 12-bit 845–885 0.06–281 12-bit

4
C. Macfarlane et al. Remote Sensing Applications: Society and Environment 6 (2017) 1–14

2.2.3. Atmospheric correction same density of vegetation) converge on the origin (Huete, 1988). The
To atmospherically correct each scene we used the simple DOS resulting SR and NDVI scenes were also saved as 8 bit GeoTIFF files.
(Chavez, 1988) method, which assumes the existence of dark objects
with zero or small surface reflectance throughout a Landsat scene and a 2.2.5. Inter-sensor comparison of SVIs
horizontally homogeneous atmosphere (Song et al., 2001). For situa- To ensure that isolines of NDVI and SR from different Landsat
tions in which absolute reflectance is not required simple methods such sensors represented the same vegetation density, we undertook a
as DOS are recommended (Chavez, 1988; Song et al., 2001). In this comparison of NDVI and SR scenes from the different Landsat sensors.
approach, clear, dark water and deep shadows are commonly used A full comparison of all Landsat sensors is difficult owing to the limited
targets (Ahern et al., 1977; Gordon, 1978; Chavez, 1988; Chang et al., availability of scenes prior to 1980 and the need to mosaic early scenes
2008). We obtained the dark point from approximately 4000 pixels from multiple MSS sensors in some cases. For inter-sensor compar-
(30 m resolution) located over deep oligotrophic water within four isons, scenes from the study series were divided into four groups
long-established ( > 40 years) reservoirs (see Fig. 1). We did not set the following Chander et al. (2009) and Flood (2014):
dark point to zero as in some studies (e.g. Chavez, 1988) in recognition
of the fact that the reflectance of even very dark objects is not exactly MSS: Landsat Multispectral Scanner 1–5.
zero (Chavez, 1989). The surface reflectance for deep water is less in TM: Landsat Thematic Mapper 4/5.
the NIR band than the red band; typical values are 1–2% in the NIR ETM+: Landsat Enhanced Thematic Mapper 7.
band and 2–3% in the red band (Ahern et al., 1977; Groeneveld and OLI: Landsat Operational Land Imager 8.
Barz, 2010; Oke, 1987). We ‘anchored’ the dark point for all scenes at
TOA reflectance values of 1.7% and 2.9%. This served to fix the In the case of MSS and TM, scenes from the same dates were
intercept of the soil line very close to the origin of a plot of NIR versus compared (MSS4/5 and TM4/5) while for TM, ETM+ and OLI it was
red of the atmospherically corrected scenes, as required for calculating necessary to compare scenes eight days apart. We assumed that this
ratio SVIs. Moran et al. (1992) found that varying the assumed time period was small enough so that any changes on the ground would
reflectance of dark objects from 1% to 2% had little impact on errors be negligible. The scenes used for the inter-sensor comparisons are
in image-based atmospheric corrections. The adjustments to the red listed in Table S2. A total of 2260 points were sampled on a regular grid
and NIR bands are summarised in Eqs. (1) and (2): across the region avoiding data that were missing in ETM+ scenes
Δρ N = (ρ N *− ρ N) owing to failure of the Scan Line Corrector (Chander et al., 2009) from
(1)
May 2003 onwards. Robust least squares linear regression (Holland
Δρ R = (ρ R *− ρ R) (2) and Welsch, 1977) was used to construct relationships between ETM+
and OLI, between TM and ETM+ and between MSS and TM. To ensure
where ΔρN and ΔρR are the adjustments required in the NIR and red that SR and NDVI isolines from the different Landsat sensors were
bands in each scene, ρN*(1.7%) and ρR*(2.9%) are the target values for comparable, these regressions were applied to all NDVI and SR scenes
reflectance of deep water, and ρN and ρR are the original values for to create scenes in the following order: the ETM+/OLI regression was
reflectance for deep water in each scene. applied to all (non-OLI) scenes, the TM/ETM+ regression was then
also applied to all MSS and TM scenes, and the MSS/TM regression
2.2.4. Adjustment of the soil line and calculation of SVIs was then also applied only to MSS scenes. Where these regression lines
Following atmospheric correction, we adjusted the NIR band of all passed close to the origin (‘0,0’ in the case of NDVI regressions or ‘1,1’
scenes so that, on average, bare soil pixels had a SR of one and a NDVI in the case of SR regressions), the regressions were forced through the
of zero. The soil line was obtained from known locations of bare ground origin on the basis that the atmospheric correction and soil line
based on a > 40 year record of bauxite mining operations in the study adjustments already applied should fix the intercept at the origin in a
region (Alcoa, unpublished data). Unlike the selection of bare ground comparison of SVIs from different sensors. Adjustment of the dark
based on features within the imagery itself, which can result in sparsely point and soil line should result in bare soil having a SR=1 and
or variably vegetated regions being included in the soil line, this long NDVI=0 for all scenes regardless of the sensor.
and spatially-explicit record allowed the selection of precisely known
locations that had been cleared of vegetation and excavated but not yet 2.3. Ground-based estimation of LAI and calibration of vegetation
revegetated, thereby reliably detecting only 100% bare ground to indices
calculate the soil line. The excavated lateritic soil profile, along with
cleared but unexcavated areas with intact topsoil, provided a range of We compiled a dataset of 1228 indirect ground-based plot mea-
digital numbers in the red and NIR bands that were typical of the study surements of overstorey and understorey photosynthetically active (i.e.
area. Further controls were employed to eliminate mixed pixels from green vegetation) LAI in the northern jarrah forest collected in the
the perimeters of mine pits and from errors arising from mining period 1991–2015 (Alcoa, unpublished data; Wallace, 1992; Wallace,
records or registration uncertainties. The slope of the soil line was 1996; Macfarlane et al., 2007a; Macfarlane et al., 2007b; Boer et al.,
estimated using robust least squares linear regression (MATLAB's 2008; McCaw, 2011; Mauger et al., 2013). The dataset comprised
‘robustfit’ function with default settings forced through the origin) measurements using either hemispherical photograph (HP, Rich,
based on iteratively reweighted least-squares (Holland and Welsch, 1990) or digital cover photography (DCP, Macfarlane et al., 2007b).
1977), and the reflectances of the NIR band adjusted by dividing by the These indirect methods have been calibrated against direct measure-
slope of the calculated soil line. ments of LAI in several studies (Carbon et al., 1979; Whitford et al.,
Following atmospheric correction and adjustment of soil line to the 1995; Macfarlane et al., 2007a, 2007b). Understorey LAI for these plots
1:1 line in NIR-red space, SR (Jordan, 1969) and NDVI (Rouse et al., was estimated using a variety of methods. Direct field-based observa-
1974) were calculated as follows from the DNs of the NIR and red tions of understorey LAI are rare (Macfarlane et al., 2010) but
bands (Eqs. (3) and (4)): estimates of understorey above-ground biomass (BAG, ton ha−1),
fractional foliage cover (ff, 0 to 1) or cumulative foliage cover (cf, %)
SR=ρ N / ρ R (3)
are more often available. We calculated understorey LAI (LUS) from
NDVI =(ρ N – ρ R)/(ρ N + ρ R) (4) these data using the following equations:
LUS = 0.193BAG (5)
Both are ratio indices that assume the slope and intercept of the soil
line are one and zero and that vegetation isolines (representing the (Macfarlane et al., 2010)

5
C. Macfarlane et al. Remote Sensing Applications: Society and Environment 6 (2017) 1–14

LUS = f f (1 − Φ ) ln(Φ )/ k (6) drop interaction terms and main effects with the largest P-values (non-
significant) until only significant (P < 0.05) terms remained in the
(Macfarlane et al., 2012)
model (Sokal and Rohlf, 1995).
ln (BAG ) = 0.94ln (cf ) − 2.73 (7)
3. Results
(Grigg and Corless, 2009) where Φ is crown porosity (assumed to be
0.3 for understorey) and k is the light extinction coefficient (assumed to 3.1. Atmospheric correction and soil line adjustment of Landsat
be 0.5). It was necessary to use both Eqs. (7) and (5) to derive Lus from scenes
cf. Cumulative foliage cover is the summed fractional foliage cover of
each species individual cover within a plot, i.e. no allowance is made for Red and NIR TOA reflectances of the dark pixels were very similar
overlapping of plant canopies. for scenes from the same sensor (Fig. 3) with the exception of Landsat
Records within the dataset were removed if they were obtained 1–3 MSS scenes, which were highly variable. The TOA reflectances of
from plots smaller than 0.16 ha, were derived from less than five HP dark pixels were also similar across some sensors, especially Landsat 4/
images or 16 DCP images per plot, were located close to abrupt 5 MSS/TM and Landsat 7 ETM+. Similarly, within a sensor the slope of
boundaries in forest cover such as clearings for roads, exact (i.e., the soil line was very similar across scenes (Fig. 3). However, the
+/−5 m) GPS co-ordinates were uncertain, or where any accompanying characteristic slope of the soil line for different sensors varied from 0.9
information on understorey biomass, cover or leaf area was unavail- to 1.2, a difference of up to 25% between sensors. The OLI sensor, from
able. A final dataset of 234 measurements was thereby attained, which the reference scene was taken, had similar slopes to the TM and
spanning the period 1996–2015. The annual time series of Landsat ETM+ sensors, while the slopes from the Landsat 4/5 MSS scenes were
imagery targeted summer and most ground observations of LAI were much smaller and those from the Landsat 1–3 MSS scenes were often
obtained in summer (150 observations), but some ground measure- larger.. The 1983 scene (day of year 61) had an unusually small soil line
ments of LAI were obtained in late spring (32 observations) or early slope, which was also observed in another scene tested for that year
autumn (52 observations; none were obtained in winter). We obtained (1983, day of year 29).
additional scenes (see Table S3) for calibration purposes, which were As expected from the similar spectral characteristics of sensors
within two months of the ground measurement's acquisition date. (Chander et al., 2009) there was little difference in NDVI or SR
These were processed identically to scenes from the main time-series. between the TM and ETM+ sensors for scenes only 8 days apart
The 234 ground measurements were calibrated against Landsat scenes (Figs. 4 and 5; Tables 3, 4); there was more difference between NDVI
from 17 different dates. and SR from ETM+ and OLI and the relationship was slightly weaker
(R2=0.92) than that between TM and ETM+ (R2=0.96). There were
2.4. Catchment LAI and streamflow large differences between the TM and MSS sensors although the
correlation between the two sensors was still high (R2=0.86–0.87). In
Annual streamflow and rainfall data were collated from the Western the case of the comparison of SR values this was reflected in a large
Australian Department of Water summary records (http://wir.water. slope (1.66) although the regression line passed close to the origin. In
wa.gov.au/; last accessed 19/Feb/2016) for a total of 30 gauged the case of NDVI the slope was closer to one (1.08) but there was a
catchments (Fig. 1, Table 2). Rainfall was typically measured at the large intercept (0.088). The relationship between ETM+ and OLI for
stream gauge, however, for larger catchments or for catchments with NDVI also had a large intercept, whereas the relationship between TM
known strong rainfall gradients (Dingo Rd, Hillview Farm, McKnoes, and ETM+ for NDVI passed very close to the origin.
North Rd, Sandalwood, Vardi Rd), either a more centrally located
gauge or an average of several gauges within the catchment was used. 3.2. LAI regression model
Only catchments with at least nine years of streamflow and rainfall
records in the period 1973–2014 were considered. All 30 catchments Given the large intercepts on the NDVI relationship between TM
were from higher rainfall parts of the region ( > approx. 1000 mm/year and MSS data, and between OLI and ETM+ data, we chose to calibrate
annual average) where groundwater connection to streams in the ground-based LAI measurements against SR values rather than NDVI
catchment has historically been strong. Some catchments (Bennett, for the purpose of constructing a long time series of LAI. However, we
Lewis, Del Park, Jack Rocks, More Seldon Seen, Warren) experienced present calibrations of LAI against NDVI for comparison, noting that
disconnection in groundwater following a record drought year in 2010, no MSS scenes were used in the calibration. Total LAI (overstorey and
causing a fundamental change in flow response in these catchments understorey combined) in the ground-based measurements ranged
(Kinal and Stoneman, 2012), and subsequent flow years were not from almost zero to nearly four and the accompanying SR ranged
included in the analysis (Table 2). The catchments encompass forest between 1 and 7 (NDVI ranged between 0.1 and 0.8). A robust
areas with varying degrees of disturbance in the form of timber regression of total LAI against SR, although linear (Fig. 6), was
harvesting, fuel reduction burning, bauxite mining and restoration, relatively poor (R2=0.59; n=234, RMSE=0.74, p < 0.001). The regres-
or experimental forest thinning including one catchment that was sion of total LAI against NDVI (Fig. 6) yielded a stronger relationship
completely cleared to pasture (Wights). To identify whether a mining (R2=0.69, n=234, RMSE=0.50, p < 0.001) that was best described by a
land use was associated with a significantly different flow response, power function. The intercept did not differ significantly from zero and
each year of flow record in the catchments dataset (n=861) was was not fitted in the final model; however, despite the stronger
assigned to either a mined or unmined category, with the first year of correlation, the relationship of LAI to NDVI appeared to overestimate
mining deemed to be one year after initial clearing. Once mining had very small values of LAI – far more so than the regression of LAI
occurred, the catchment was classed as mined for the remainder of the against SR. In the ground-based LAI plots, overstorey accounted for
available record. Two catchments with mining disturbance on the more than two-thirds of total plot LAI in over 85% of plots (see Table
catchment divide and constituting only 2% of the catchment area were S4).
classed as unmined (Table 2).
We used generalized linear models implemented in Minitab 17 3.3. Catchment LAI and streamflow
(Minitab Inc., Pennsylvania, USA) to test for the response of annual
streamflow to annual rainfall, annual LAI spatially averaged for the Changes in catchment average LAI since 1973 in the series
whole catchment, and interactions. Mining was included as a random developed in this study, and differences in LAI between catchments,
factor. Subsequently, backwards elimination was used to sequentially are consistent with a range of treatments applied and measured forest

6
C. Macfarlane et al. Remote Sensing Applications: Society and Environment 6 (2017) 1–14

Table 2
Characteristics of catchments used in streamflow prediction.

Catchment Area (km2) Record start Record end Record length % mineda Start of mining End of miningb Streamflowc Rainfallc (mm)
(years) (mm)

31 Mile Bk 11 1986 2014 22 2 – – 109 984

Bates 2.3 1989 2014 26 0 – – 131 1095
Bussell Brook 38 1978 1998 21 0 – – 107 1010
Davis Brook 66 1973 2001 29 0 – – 65 1050
Dingo Rd 147 1973 2014 42 0 – – 198 1149
Higgens 0.7 1978 1998 21 0 – – 86 1191
Hillview Farm 17 1977 2014 38 0 – – 213 1189
Jack Rocks 55 1982 2010 21 2 – – 86 1000
Jones 0.7 1978 1998 21 0 – – 52 1161
Salmon 0.8 1975 2013 29 0 – – 105 1092
Sandalwood 116 1981 2014 34 0 – – 196 1131
Waterfall Gully 9 1975 2014 40 0 – – 195 1074
Wights 0.9 1974 2006 33 0 – – 372 1109
Worsley 32 1973 2014 42 0 – – 151 1014
Banksiadale 3.7 2002 2014 13 56 1982 1993 91 1171
Bennetts 0.7 1978 2010 33 48 1989 1992 242 1144
Conjuranup 39 1973 1991 14 43 1976 2007 206 1219
Del Park 1.3 1975 2010 36 37 1975 1989 156 1148
Hansens 0.7 1978 2014 30 62 2000 2006 168 1161
Lewis 1.9 1978 2010 33 49 1997 2003 108 1144
Marrinup Brook 46 1973 2014 42 16 1971 1992 186 1175
McKnoes 24 1980 2014 26 32 1980 2001 312 1202
More Seldom Seen 3.3 1975 2010 36 57 1970 1994 181 1076
Mt William 4.4 1986 2014 20 22 1990 2000 334 1252
North Rd 80 1984 2014 24 32 1992 2013 71 1088
Seldom Seen 7.4 1975 2014 40 33 1968 2001 174 1048
Skeleton Rd 16 1989 1997 9 30 1991 2007 161 1125
Vardi Rd 81 1982 2014 33 22 1968 2000 111 1028
Warren 0.8 1978 2010 33 40 1989 1992 164 1137
Waterous 14 1975 2014 40 7 1994 2000 262 1188

a
– Maximum area mined in the period of flow record.
b
– Mining has ended when all restoration is complete.
c
– Mean for period of flow record.

cover reported previously for six of the catchments included in this green vegetation and 3) adjusted the SR values of earlier Landsat
study (Fig. 7). We therefore have confidence that both the spatial and sensors to ensure that ratio index isolines were consistent with the
temporal estimates of LAI developed here are realistic. most recent imagery from the Landsat 8 OLI. The method was applied
Both LAI and annual rainfall had a highly significant effect on to Landsat scenes that were converted to TOA reflectance but could
streamflow (GLM, p < 0.001) in the model developed for all 30 equally be applied to any time series of imagery, even starting with raw
catchments, with an inverse relationship between streamflow and digital numbers. As a result the method could have application to
LAI and a positive relationship between streamflow and rainfall integrating time series from sensors other than Landsat and for which
(Fig. 8). A highly significant LAI-rainfall interaction was also evident, calibrations to reflectance or radiance are not available.
however, mining (excavation) was not a significant factor (p=0.81). The A significant advantage of our approach is that it used features of
final model predicting streamflow was given by Eq. (8): the landscape for relative radiometric calibration that represented the
environment for which the time series was required. The landscape
Q = −434 + 111 L + 0.755 P –0.226 L (n = 861, R2 = 0.58) (8)
features that we used to standardize the time series, surface water
where Q is annual streamflow (mm), L is average catchment LAI for the reservoirs and bare ground, also represented the landscape processes
series year scene and P is annual rainfall (mm). that we wished to investigate: the study aimed to quantify the
relationship between clearing and excavation due to mining and
4. Discussion surface water run-off. Manual selection of invariant targets is labour
intensive and the spectral characteristics of targets may change over
The availability of much improved MSS Landsat imagery from the time, either if the land use changes or if targets are not truly bare such
USGS (Woodcock et al., 2008; Wulder et al., 2012) since 2010, with that sparse vegetation may appear (Hall et al., 1991; Schott et al.,
improved radiometric calibrations, improved geometric corrections 1988). The long-term GIS record of clearing, mining and revegetation
(Markham and Helder, 2012) and greater consistency with later made it possible to automate the identification of ‘true’ bare ground
Landsat sensors, has renewed interest in the integration of MSS targets extending over 40 years and to construct bare soil lines for each
imagery into long time series for the purposes of studying forest scene. Used together with deep, oligotrophic water sources in surface-
dynamics (e.g. Gómez et al., 2011; Pflugmacher et al., 2012). Unlike water reservoirs we were able to standardize Landsat scenes since they
other recent studies (Powell et al., 2008; Gómez et al., 2011; first became available in 1972. Such a detailed record of land clearing
Pflugmacher et al., 2012), we successfully produced a long annual time may not be available in many cases, and in arid regions permanent
series of LAI for the northern jarrah forest which integrated Landsat water bodies might not be available as targets for dark point adjust-
MSS sensor imagery with later TM, ETM+ and OLI imagery, using a ment. Identifying bare ground targets retrospectively involves great
simple ratio rather than a tasselled cap approach. Our simple approach challenges if precise land use records are not available, hence the
1) adjusted the dark point (atmospheric correction) to ensure that ratio development of automated statistical methods to identify pseudo-
index isolines converge at the origin (Donohue et al., 2008), 2) adjusted invariant targets (e.g. Canty et al., 2004; Paolini et al., 2006; Schmidt
the soil line to ensure that the SR=1 isoline represented pixels with no et al., 2008).

7
C. Macfarlane et al. Remote Sensing Applications: Society and Environment 6 (2017) 1–14

Fig. 3. Dark points (% TOA reflectance) and slope of the soil line from 1973 to 2015 for each Landsat sensor (MSS=Multispectral Scanner 1–4; TM=Thematic Mapper 5; ETM
+=Enhanced Thematic Mapper 7; OLI=Operational Land Imager 8).

Within the TM, ETM+ and OLI sensor groups the slopes of bare soil istics of the two sensors (Table 1) despite the eight day difference in
lines were similar, although there were consistent differences between acquisition date. Numerous comparisons of MSS and TM have been
sensors (Fig. 3). The extreme consistency in the bare soil lines from the done in the past (e.g. Price, 1987; Royer et al., 1987; Huete et al., 2002;
long (i.e. 1988–2011) Landsat 5 TM record gives us great confidence in Steven et al., 2003; Markham and Helder, 2012;) and the similarity of
our automated selection of bare soil sites. Together with atmospheric TM and ETM+ has been previously documented (e.g. Tiellet et al.,
effects, it is likely that the spectral characteristics of the different 2001; Steven et al., 2003), and confirmed in our study. There was also
sensors (Table 1) is contributing to the observed differences in the dark good correspondence between ETM+ and OLI derived values of SVIs,
points of scenes from different Landsat sensors, especially between again despite the eight day difference in acquisition date. Band 5 (NIR)
MSS scenes and later sensors. Spectral differences between the sensor of the OLI sensor is much narrower than the corresponding band of TM
groups were most evident in the red band (Fig. 3), despite there being or ETM+ to avoid an atmospheric water absorption feature at 0.825 µm
particularly poor correspondence of TM Band 4 (NIR) with either (Irons et al., 2012) and recent studies have noted small differences in
MSS1-3 bands 6 or 7 or MSS4-5 bands 3 or 4. We selected band 3 from NDVI calculated from ETM+ versus OLI as a result. Li et al. (2014)
MSS4-5 to represent NIR (band 6 from MSS1-3) based on the more found that NDVI from OLI was generally larger than that from ETM+
certain cross-calibration of this band compared to the longer wave- by about 0.05. Similarly, Flood (2014) found that, prior to any
length band (Markham and Helder, 2012). Galvao et al. (1999) found adjustments to reflectance of individual bands, OLI yielded values of
that variations in the width of the red band had little impact on NDVI that were 5% larger than ETM+. Xu and Guo (2014) concluded,
calculated NDVI provided that the band did not extend beyond based on simulations, that NDVI from OLI would be marginally larger
0.69 µm. The variation of dark point reflectances and bare soil slopes than that from ETM+ but that this would be most evident at small
between the various MSS sensors may relate to this uncertainty values of NDVI. Our study corroborates these previous studies in that
surrounding the calibration of the MSS sensors: Hall et al. (1991) values of NDVI from OLI were generally larger than those from ETM+
reported that, post-launch, calibrations of Landsat 2–5 MSS only (Table 3 and Fig. 4), and the difference was slightly larger at small
agreed to within 12% of one another and tended to diverge further values. There was clearly more variation (and a correspondingly
with age. Recent cross-calibration (Chander et al., 2009) has reduced smaller R2 value) in the relationship between SVIs from OLI and
uncertainty in the three shorter bands to ~4%, while the longest ETM+, than between SVIs from ETM+ and TM, which presumably
wavelength band still has much greater uncertainty. results from greater atmospheric interference in the NIR band of ETM
The effect of spectral variations between different Landsat sensors + and TM compared to OLI.
was also evident in the inter-sensor comparison of SR and NDVI values However, the greatest differences were observed between SVIs
(Figs. 4 and 5). There was a strong relationship between TM and ETM+ derived from MSS and TM. This was despite the data being obtained
derived values of SVIs, reflecting the very similar spectral character- from the same Landsat satellite on the same date and time. The slope of

8
C. Macfarlane et al. Remote Sensing Applications: Society and Environment 6 (2017) 1–14

Fig. 4. Density plots of NDVI comparing different sensor groups (MSS=Landsat

Multispectral Scanner 5; TM=Landsat Thematic Mapper 5; ETM+=Landsat Enhanced
Thematic Mapper 7; OLI=Landsat Operational Land Imager 8). Colours represent the
frequency of X-Y combinations; any X-Y combinations that occur more than 16 times are Fig. 5. Density plots of simple ratio (SR) comparing different sensor groups
shown as yellow. (For interpretation of the references to color in this figure legend, the (MSS=Landsat Multispectral Scanner 5; TM=Landsat Thematic Mapper 5; ETM
reader is referred to the web version of this article). +=Landsat Enhanced Thematic Mapper 7; OLI=Landsat Operational Land Imager 8).
Colours represent the frequency of X-Y combinations; any X-Y combinations that occur
the regression for SR was much larger than one (Table 4) but the more than 16 times are shown as yellow. (For interpretation of the references to color in
regression line passed close to the origin; it is important for the this figure legend, the reader is referred to the web version of this article).

regression to pass through the origin because the adjustment of the

dark point and soil line should result in bare soil having a SR=1 and homogeneous forested area, there were areas of annual pasture that
NDVI=0 for all scenes regardless of the sensor. The situation was much usually had small NDVI and SR values during the dry summer period
worse for NDVI: the regression line clearly did not pass close to the that we targeted. In Fig. 4, these areas are represented by the small
origin and all values of NDVI from TM were larger than those from the cluster of points with NDVI between 0.1 and 0.2. The fact that SR
corresponding MSS scene. Hence, the SR may offer great advantages allowed us to create a single translation equation across the land types
over NDVI for integrating MSS imagery into LTS using ratio indices suggests that SR may offer distinct advantages over NDVI for integrat-
across different land classes. Although our study targeted a relatively ing different sensors to create long time series. It has previously been

9
C. Macfarlane et al. Remote Sensing Applications: Society and Environment 6 (2017) 1–14

Table 3 observed that separate equations may be needed for distinct land
Regression parameters from comparisons of NDVI from four scenes from each sensor classes (Yoshioka, 2003) owing to the interactions between the spectral
group after standardizing for dark point and bare soil line.
characteristics of surface vegetation and soil components, and sensor-
Y X Slope p (slope) Intercept p (int) R2 D.F. specific spectral band characteristics (van Leeuwen et al., 2006).
In our study the relationship between LAI and NDVI was best
TM MSS 1.08 < 0.001 0.088 < 0.001 0.87 9038 described by a non-linear relationship while the relationship between
ETM+ TM 1.02 < 0.001 0.013 < 0.001 0.97 9038
LAI and SR was best described by a linear relationship. A strength of
OLI ETM+ 0.97 < 0.001 0.048 < 0.001 0.95 9038
our study is the large number of ground observations of LAI, and the
calibration of SVIs against total LAI – not just overstorey LAI. The wide
Table 4 variation in sensors used, types and densities of vegetation studied,
Regression parameters from comparisons of ‘SR – 1’ from four scenes from each sensor ground methods employed, definitions of LAI applied, and pre-
group after standardizing for dark point and bare soil line. processing of imagery makes a comparison of LAI-SVI relationships
with other studies difficult. Our results agree with previous studies that
Y X Slope p (slope) R2 D.F.
found SR was linearly correlated with LAI (McVicar et al., 1996; Coops
TM MSS 1.66 < 0.001 0.86 9038 et al., 1997; Lu et al., 2003; Zarate-Valdez et al., 2012). Both linear and
ETM+ TM 1.10 < 0.001 0.96 9038 non-linear relationships have been reported between LAI and NDVI
OLI ETM+ 1.15 < 0.001 0.94 9038 (Asrar et al., 1985; Gamon et al., 1995; McVicar et al., 1996; Coops
et al., 1997; Qi et al., 2000; Lu et al., 2003; Zarate-Valdez et al., 2012;
Tripathi et al., 2014). Carlson and Ripley (1997) argued that the
relationship between NDVI and LAI would be nearly linear at low LAI
values, where substantial bare ground may be visible from above, but
become increasingly curved as fractional vegetation cover approached
100%. The coefficient of determination (0.69) we obtained for the
regression of LAI versus NDVI is less than that obtained by Boeret al.
(2008; 0.87) for similar vegetation using Landsat TM and similar
ground methods, and the coefficient of determination for our regres-
sion of LAI versus SR was also relatively small (0.59). However, the
relationships obtained by Boer et al. (2008) were based on a single
Landsat scene, and was restricted to a smaller area of the northern
jarrah forest and a single season of ground-based data. That study also
only considered overstorey measurements of LAI and ignored under-
storey. Whilst our correlation coefficients were smaller, the large
number of ground observations and temporal span of observation
across multiple, cross-calibrated sensors makes our calibration very
robust for application across the entire LTS. Despite the smaller
correlation coefficient of our SR-based regression, we used SR to
model LAI from the LTS owing to the large positive intercept on the
regressions comparing NDVI between MSS and TM, and between ETM
+ and OLI, which would have resulted in NDVI isolines from the
different sensors not representing the same vegetation amount.
There are many sources of uncertainty in ground-based estimates of
LAI (Richardson et al., 2011). There are few accurate estimates of the
LAI of understorey available and we estimated this from a variety of
other measurements based on numerous assumptions. This is likely to
be a significant source of uncertainty, even greater than the uncertainty
in estimates of overstorey LAI, which are likely to exceed 10% and
could be as large as 30% (Chen et al., 1997; Raupach et al., 2005).
Overstorey estimates of LAI from DCP require an assumed value of the
light extinction coefficient, which assumes a fixed leaf angle distribu-
tion. However, leaf angles can vary with wind speed and strong winds
can increase measured cover in erectophile canopies by up to 10%
(Macfarlane and Landman, 2002). Taking account of wind speed (e.g.
McVicar et al., 2008) at the time of both ground measurement and
scene acquisition could help to refine calibrations of remotely-sensed
indices against ground-based observations. Scenes were also some-
times obtained up to two months apart from ground observations
which could lead to small variations in the LAI between sampling
times. There will necessarily be greater uncertainty associated with
estimates of LAI from earlier MSS sensors owing to the coarser spatial
resolution and broader bands (Table 1), and also because the record of
ground calibration data did not extend that far back in time.
Fig. 6. Relationship between NDVI and SR from Landsat imagery and total stand LAI
The model of annual streamflow developed in this study is one of
from ground measurements for the study region.
few empirical attempts to provide a single framework for under-
standing the relationship between forest density and streamflow, as
recommended by Robinson et al. (1997) almost two decades ago. While
the strength and direction of the effects of rainfall and LAI on

10
C. Macfarlane et al. Remote Sensing Applications: Society and Environment 6 (2017) 1–14

Fig. 7. Changes in catchment average LAI over time under a range of applied treatments: Hansens (a) was heavily thinned in 1985/86, reducing canopy cover from 50 to 60–14% (basal
area 27–7 m2/ha) and allowed to regenerate (Ruprecht et al., 1991) before being mined for bauxite from 2002 and subsequently rehabilitated. Higgins (b) and Jones (c) were subject to
experimental thinning in 1988/89 reducing forest basal area from approximately 40 m2/ha to 15 and 17 m2/ha, respectively (Robinson et al., 1997). Both catchments were mined for
bauxite from 2002 and 2001, respectively, and rehabilitated. A fully forested catchment, Wights (d), was completely cleared for pasture development in 1977 (Peck and Williamson,
1987). More Seldom Seen (e) was subject to bauxite mining and rehabilitation over the period 1970–1994 with a peak in area cleared but not rehabilitated in 1982 (Bari and Ruprecht,
2003). A control catchment Waterfall Gully (f) was severely affected by the soil pathogen Phytophthora cinnamomi, which causes death in jarrah and a range of forest species, reducing
forest cover (Batini et al., 1980).

streamflow are not unexpected, the absence of a significant mining- storage as a catchment's ‘memory. Without specific representation of
specific response suggests consistency across different land use types in the groundwater storage component, our model was therefore likely to
the study region, which simplifies parameterization of catchment be unable to fully account for such longer-term trends leading to the
runoff or process-based models. Streamflows in the region exhibit relatively low fraction of variation explained by the GLM.
variable rainfall-runoff responses described as a type of ‘non-statio- The model is consistent with the eco-hydrological equilibrium
narity’ (Petrone et al., 2010; Hughes and Vaze, 2015) linked to longer- theory (Eagleson, 1982) whereby vegetation at equilibrium with
term changes in catchment groundwater storage (Hughes et al., 2012). climate evapotranspires almost all rainfall leaving minimal excess for
Sustained reductions in groundwater storage have been observed after other components of the water balance including runoff. Schofield et al.
drought years, leading Hughes et al. (2012) to describe groundwater (1989), for example, observed that dense jarrah forest growing on more

11
C. Macfarlane et al. Remote Sensing Applications: Society and Environment 6 (2017) 1–14

hydrological prediction in the region by incorporating long term

changes in catchment forest LAI.

Acknowledgements

AHG and MID are paid employees of Alcoa of Australia Ltd. CM

was partly supported by the Australian Supersite Network, part of the
Terrestrial Ecosystem Research Network (www.tern.org.au). We thank
Tim McVicar for comments on an early version of this manuscript that
greatly improved the manuscript.

Appendix A. Supplementary material

Supplementary data associated with this article can be found in the

online version at doi:10.1016/j.rsase.2017.01.006.

References

Ahern, F.J., Goodenough, D.G., Jain, S.C. Rac, V.R., Rochon, G., 1977. Use of clear lakes
as standard reflectors for atmospheric measurements. In: Proceedings of the 11th
International Symposium on Remote Sensing of Environment, Ann Arbor, MI, pp.
731–775.
Aitken, S.N., Yeaman, S., Holliday, J.A., Wang, T., Curtis-McLane, S., 2008. Adaptation,
migration or extirpation: climate change outcomes for tree populations. Evolut. Appl.
1, 95–111.
Asrar, G., Kanemasu, E.T., Yoshida, M., 1985. Estimates of leaf area index from spectral
reflectance of wheat under different cultural practices and solar angle. Remote Sens.
Environ. 17, 1–11.
Baldocchi, D.D., Harley, P.C., 1995. Scaling carbon dioxide and water vapour exchange
from leaf to canopy in a deciduous forest. II. Model testing and application. Plant
Cell Environ. 18, 1157–1173.
Baret, F., Guyot, G., 1991. Potentials and limits of vegetation indices for LAI and APAR
assessment. Remote Sens. Environ. 35, 161–173.
Baret, F., Jacquemoud, S., Hancoq, J.F., 1993. About the soil line concept in remote
sensing. Adv. Space Res. 13, 281–284.
Bari, M.A., Ruprecht, J.K., 2003. Water yield response to land use change in south–west
Western Australia, Department of Environment, Salinity and Land Use Impacts
Fig. 8. Annual streamflow in relation to (a) catchment average LAI for that year and (b) Series Report no. SLUI 31.
annual rainfall, for 30 catchments in the northern jarrah forest of south-western Bates, B.C., Hope, P., Ryan, B., Smith, I., Charles, S., 2008. Key findings from the Indian
Australia (Fig. 1; Table 2). Ocean Climate Initiative and their impact on policy development in Australia. Clim.
Change 89, 339–354.
Batini, F.E., Black, R.E., Byrne, J., Clifford, P.J., 1980. An examination of the effects of
favourable sites had the capacity to evapotranspire virtually all annual changes in catchment condition on water yield in the Wungong Catchment, Western
precipitation. For a long-term annual average rainfall of 1200 mm, our Australia. Aust. For. Res. 10, 29–38.
Bell, D.T., Heddle, E.M., 1989. Floristic, morphologic and vegetational diversity. In: Dell,
model predicts zero streamflow when catchment average LAI reaches B., Havel, J.J., Malajczuk, N. (Eds.), The Jarrah Forest: A Complex Mediterranean
2.9. This matches reasonably well with a predicted equilibrium LAI of Ecosystem. Kluwer Academic Publishers, Dordrecht, 53–66.
3.2 for the same rainfall using the relationship developed by Ellis and Boer, M.M., Macfarlane, C., Norris, J., Sadler, R.J., Wallace, J., Grierson, P., 2008.
Mapping burned areas and burn severity patterns in SW Australian eucalypt forest
Hatton (2008). Together with the close match between changes in using remotely-sensed changes in leaf area index. Remote Sens. Environ. 112,
catchment LAI in our LTS and changes in catchment forest cover due to 4358–4369.
experimental treatments reported by others (Fig. 7), we therefore Bréda, N.J.J., 2003. Ground-based measurements of leaf area index: a review of
methods, instruments and current controversies. J. Exp. Bot. 54, 2403–2417.
conclude that the long time-series of LAI developed in this study is Canadell, J.G., Raupach, M.R., 2008. Managing forests for climate change mitigation.
robust and potentially valuable for a range of hydrological applications, Science 320, 1456–1457.
extending and enhancing recent investigations of streamflow in the Canty, M.J., Nielson, A.A., Schmidt, M., 2004. Automatic radiometric normalization of
multitemporal satellite imagery. Remote Sens. Environ. 91, 441–451.
region that have focused on changes in rainfall only (Petrone et al.,
Carbon, B.A., Bartle, G.A., Murray, A.M., 1979. Leaf area index of some eucalypt forests
2010; Hughes et al., 2012; Silberstein et al., 2012; but see Smettem in south-west Australia. Aust. For. Res. 9, 323–326.
et al., 2013). Carlson, T.N., Ripley, D.A., 1997. On the relation between NDVI, fractional vegetation
cover, and leaf area index. Remote Sens. Environ. 62, 241–252.
Chander, G., Markham, B.L., Helder, D.L., 2009. Summary of current radiometric
calibration coefficients for Landsat MSS, TM, ETM+, and EO-1 ALI sensors. Remote
5. Conclusion Sens. Environ. 113, 893–903.
Chang, J., Clay, D.E., Leigh, L., Aaron, D., Dalsted, K., Volz, M., 2008. Evaluating
In this study, we developed a standardized annual time-series of modified atmospheric correction methods for Landsat imagery: image-based and
model-based calibration methods. Commun. Soil Sci. Plant Anal. 39, 1532–1545.
LAI for the northern jarrah forest of south-western Australia using the Chavez, P.S., 1988. An improved dark-object subtraction technique for atmospheric
44 year time series of level 1 terrain corrected Landsat imagery sourced scattering correction of multispectral data. Remote Sens. Environ. 24, 459–479.
from the USGS. A feature of the series was the successful integration of Chavez, P.S., 1989. Radiometric calibration of Landsat Thematic Mapper multispectral
images. Photogramm. Eng. Remote Sens. 55, 1285–1294.
Landsat MSS sensor imagery with later TM, ETM+ and OLI imagery,
Chen, J.M., Rich, P.M., Gower, S.T., Norman, J.M., Plummer, S., 1997. Leaf area index of
using a simple ratio rather than a tasselled cap approach. Our generic boreal forests: theory, techniques and measurements. J. Geophys. Res. 102 (D24),
approach was based on atmospheric correction using a dark point, 29429–29443.
Cohen, W.B., Goward, S.N., 2004. Landsat's role in ecological applications of remote
adjustment of the bare soil line and inter-sensor calibration of
sensing. Bioscience 6, 535–545.
vegetation index isolines, and could be applied to other sources of Coops, N., Delahaye, A., Pook, E., 1997. Estimation of eucalypt forest leaf area index on
imagery. The time-series of LAI was applied to develop a simple the south coast of New South Wales using Landsat MSS data. Aust. J. Bot. 45,
relationship between rainfall, forest LAI as affected by land use 757–769.
Donohue, R.J., Roderick, M.L., McVicar, T.R., 2008. Deriving consistent long-term
activities, and streamflow. The time-series has potential to improve

12
C. Macfarlane et al. Remote Sensing Applications: Society and Environment 6 (2017) 1–14

vegetation information from AVHRR reflectance data using a cover-triangle-based monitoring canopy cover at Paraburdoo and Yandicoogina and testing of
framework. Remote Sens. Environ. 112, 2938–2949. photographic and light transmission methods for assessing canopy cover. Report to
Eagleson, P.S., 1982. Ecological optimality in water-limited natural soil–vegetation Hamersley Iron Pty. Ltd. 42 pp.
systems, 1. Theory and hypothesis. Water Resour. Res. 18, 325–340. Macfarlane, C., Grigg, A., Evangelista, C., 2007a. Estimating forest leaf area using cover
Ellis, T.W., Hatton, T.J., 2008. Relating leaf area index of natural eucalypt vegetation to and fullframe fisheye photography: thinking inside the circle. Agric. For. Meteorol.
climate variables in southern Australia. Agric. Water Manag. 95, 743–747. 146, 1–12.
Flood, N., 2014. Continuity of reflectance data between Landsat-7 ETM+ and Landsat-8 Macfarlane, C., Ogden, G.N., Silberstein, R.S., 2012. Water Balance of 31 Mile Brook,
OLI, for both top-of-atmosphere and surface reflectance: a study in the Australian Western Australia. Report to the Water Foundation. CSIRO, Floreat.
landscape. Remote Sens. 6, 7952–7970. Macfarlane, C., Lardner, T., Patterson, K., Grigg, A.H., 2010. A new model for predicting
Gallo, K.P., Daughtry, C.S.T., 1987. Differences in vegetation indices for simulated understorey leaf area from biomass in eucalypt forest to test the ecohydrological
Landsat-5 MSS and TM, NOAA-9 AVHRR, and SPOT-1 sensor systems. Remote equilibrium theory. Methods Ecol. Evol. 1, 371–379.
Sens. Environ. 23, 439–452. Macfarlane, C., Hoffman, M., Eamus, D., Kerp, N., Higginson, S., McMurtrie, R., Adams,
Galvao, L.S., Vitorell, I., Filho, R.A., 1999. Effects of band positioning and bandwidth on M., 2007b. Estimation of leaf area index in eucalypt forest using digital photography.
NDVI measurements of tropical savannas. Remote Sens. Environ. 67, 181–193. Agric. For. Meteorol. 143, 176–188.
Gamon, J.A., Field, C.B., Goulden, M.L., Griffen, K.L., Hartley, A.E., Joel, G., Penuelas, Markham, B.L., Helder, D., 2012. Forty-year calibrated record of earth-surface reflected
J., Valentini, R., 1995. Relationships between NDVI, canopy structure, and radiance from Landsat: a review. Remote Sens. Environ. 122, 30–40.
photosynthesis in three Californian vegetation types. Ecol. Appl. 5, 28–41. Mauger, G.W., Grigg, A.H., Croton, J.T., 2013. Preparation of a long time series of leaf
Gardner, J.H., Bell, D.T., 2007. Bauxite mining restoration by Alcoa World Alumina area index grids for the northern jarrah forest, Western Australia. Environmental
Australia in Western Australia: social, political, historical, and environmental Department Research Bulletin 42, Alcoa of Australia.
contexts. Restor. Ecol. 15, S3–S10. McCaw, W.L., 2011. Characteristics of jarrah (Eucalyptus marginata) forest at
Gascon, F., Gastellu-Etchegorry, J.P., Lefevre-Fonollosa, M.J., Dufrene, E., 2004. Forestcheck monitoring sites in south-west Western Australia: stand structure, litter,
Retrieval of forest biophysical variables by inverting a 3-D radiative transfer model woody debris, soil and foliar nutrients. Aust. For. 74, 254–265.
and using high and very high resolution imagery. Int. J. Remote Sens. 25, McVicar, T.R., Jupp, D.L.B., Williams, N.A., 1996. Relating AVHRR Vegetation Indices
5601–5616. To Landsat TM Leaf Area Index Estimates. CSIRO Land and Water, Canberra,
Gholz, H.L., Ewel, K., Teskey, R.O., 1990. Water and forest productivity. For. Ecol. Australia, (Technical Memorandum 96/15).
Manag. 30, 1–18. McVicar, T.R., Davies, P.J., Qinke, Y., Zhang, G., 2002. An introduction to temporal-
Gómez, C., White, J.C., Wulder, M.A., 2011. Characterizing the state and processes of geographic information systems (TGIS) for assessing, monitoring and modelling
change in a dynamic forest environment using hierarchical spatio-temporal regional water and soil processes. In: McVicar, T.R., Rui, T.R., Walker, J.,
segmentation. Remote Sens. Environ. 115, 1665–1679. Fitzpatrick, R.W., Changming, L. (Eds.), Regional Water and Soil Assessment for
Gordon, H.R., 1978. Removal of atmospheric effects from satellite imagery of the ocean. Managing Sustainable Agriculture in China and Australia. ACIAR Monograph 84.
Appl. Opt. 17, 1631–1636. ACIAR, 205–223.
Grigg, A.H., Corless, E., 2009. A non-destructive method for rapid estimation of McVicar, T.R., Van Niel, T.G., Li, L.T., Roderick, M.L., Rayner, D.P., Ricciardulli, L.,
understorey biomass in bauxite mine rehabilitation. Environmental Department Donohue, R.J., 2008. Wind speed climatology and trends for Australia, 1975–2006:
Research Note 31, Alcoa of Australia. capturing the stilling phenomenon and comparison with near-surface reanalysis
Groeneveld, D.P., Barz, D.D., 2010. A robust empirical relationship for atmospheric output. Geophys. Res. Lett. 35, L20403.
scatter between red and near-infrared bands. Remote Sens. Lett. 1, 65–74. Miura, T., Huete, A., Yoshioka, H., 2006. An empirical investigation of cross-sensor
Hall, F.G., Strebel, D.E., Nickeson, J.E., Goetz, S.J., 1991. Radiometric rectification: relationships of NDVI and red/near-infrared reflectance using EO-1 Hyperion data.
toward a common radiometric response among multidate, multisensor images. Remote Sens. Environ. 100, 223–236.
Remote Sens. Environ. 35, 11–27. Moran, M.S., Jackson, R.D., Slater, P.N., Teillet, P.M., 1992. Evaluation of simplified
Havel, J.J., 1989. Land use conflicts and the emergence of multiple use. In: Dell, B., procedures for retrieval of land surface reflectance factors from satellite sensor
Havel, J.J., Malajczuk, N. (Eds.), The Jarrah Forest. Kluwer, Dordrecht, 281–314. output. Remote Sens. Environ. 14, 169–184.
Hickman, A.H., Smurthwaite, A.J., Brown, I.M., Davy, R., 1992. Bauxite mineralization Morfitt, R., Barsi, J., Levy, R., Markham, B., Micijevic, E., Ong, L., Scaramuzza, P.,
in the Darling Range, Western Australia. Report 33, Geological Survey of Western Vanderwerff, K., 2015. Landsat-8 Operational Land Imager (OLI) radiometric
Australia, Perth. performance on-orbit. Remote Sens. 7, 2208–2237.
Holland, P.W., Welsch, R.E., 1977. Robust regression using iteratively reweighted least- Oke, T.R., 1987. Boundary Layer Climates Second ed.. Routledge, London.
squares. Commun. Stat.: Theory Methods A6, 813–827. Paolini, L., Grings, F., Sobrino, J.A., Jiménez Muñoz, J.C., Karszenbaum, H., 2006.
Hostert, P., Roder, A., Hill, J., 2003. Coupling spectral unmixing and trend analysis for Radiometric correction effects in Landsat multi-date/multi-sensor change detection
monitoring of long-term vegetation dynamics in Mediterranean rangelands. Remote studies. Int. J. Remote Sens. 27, 685–704.
Sens. Environ. 87, 183–197. Peddle, D.R., Hall, F.G., LeDrew, E.F., 1999. Spectral mixture analysis and geometric
Huete, A.R., 1988. A soil adjusted vegetation index (SAVI). Remote Sens. Environ. 25, optical reflectance modeling of boreal forest biophysical structure. Remote Sens.
295–309. Environ. 67, 288–297.
Huete, A.R., Yoshioka, H., Miura, T., Kim, H.J., Gao, X., 2002. Inter-sensor calibration of Petrone, K.C., Hughes, J.D., Van Niel, T.G., Silberstein, R.P., 2010. Streamflow decline in
vegetation indices for monitoring and continuity studies of ecosystem variability. In: southwestern Australia, 1950–2008. Geophys. Res. Lett. 37, L11401.
Sobrino, J.A. (Ed.), Recent Advances in Quantitative Remote Sensing. Universiat di Pflugmacher, D., Cohen, W.B., Kennedy, R.E., 2012. Using Landsat-derived disturbance
Valancia, Spain. history (1972–2010) to predict current forest structure. Remote Sens. Environ. 122,
Hughes, J.D., Vaze, J., 2015. Non-stationarity driven by long-term change in catchment 146–165.
storage: possibilities and implications. Proc. IHAS 371, 7–12. http://dx.doi.org/ Powell, S.L., Cohen, W.B., Yang, Z., Pierce, J.D., Alberti, M., 2008. Quantification of
10.5194/piahs-371-7-2015. impervious surface in the Snohomish water resources inventory area of Western
Hughes, J.D., Petrone, K.C., Silberstein, R., 2012. Drought, groundwater storage and Washington from 1972–2006. Remote Sens. Environ. 112, 1895–1908.
stream flow decline in southwestern Australia. Geophys. Res. Lett. 39, L03408. Price, J.C., 1987. Calibration of satellite radiometers and the comparison of vegetation
Irons, J.R., Dwyer, J.L., Barsi, J.A., 2012. The next Landsat satellite: the Landsat data indices. Remote Sens. Environ. 21, 15–27.
continuity mission. Remote Sens. Environ. 122, 11–21. Qi, J., Kerr, M.S., Weltz, M., Huete, A.R., Sorooshian, S., Bryant, R., 2000. Leaf area
Jordan, C.F., 1969. Derivation of leaf-area index from quality of light on the forest floor. index estimates using remotely sensed data and BRDF models in a semiarid region.
Ecology 50, 663–666. Remote Sens. Environ. 73, 18–30.
Kauth, R.J., Thomas, G.S., 1976. The tasseled cap — A graphic description of the Raupach, M.R., Rayner, P.J., Barrett, D.J., Defries, Heimann, M., Ojima, D.S., Quegan,
spectral-temporal development of agricultural crops as seen by LANDSAT. S., Schmullius, C.C., 2005. Model–data synthesis in terrestrial carbon observation:
Symposium on Machine Processing of Remotely Sensed Data. Indiana, USA: Purdue methods, data requirements and data uncertainty specifications. Glob. Change Biol.
University of West Lafayette pp. 4B-41 to 44B-51. 11, 378–397.
Kinal, J., Stoneman, G.L., 2012. Disconnection of groundwater from surface water causes Rich, P.M., 1990. Characterizing plant canopies with hemispherical photographs.
a fundamental change in hydrology in a forested catchment in south-western Remote Sens. Rev. 5, 13–29.
Australia. J. Hydrol. 472–473, 14–24. Richardson, A.D., Bryan Dail, D., Hollinger, D.Y., 2011. Leaf area index uncertainty
Koch, J.M., 2007. Alcoa's mining and restoration process in south western Australia. estimates for model–data fusion applications. Agric. For. Meteorol. 151, 1287.
Restor. Ecol. 15, S11–S16. Robinson, J., Davies, J., van Hall, S., Bari, M., 1997. The impact of forest thinning on the
Landsberg, J.J., Waring, R.H., 1997. A generalised model of forest productivity using hydrology of three small catchments in the south west of Western Australia. Water
simplified concepts of radiation-use efficiency, carbon balance and partitioning. For. and Rivers Commission, Water Resources Technical Series No. WRT 16, Perth,
Ecol. Manag. 95, 209–228. Australia. 74p.
Li, P., Jiang, L., Feng, Z., 2014. Cross-comparison of vegetation indices derived from Rouse, J.W., Haas, R.H., Schell, J.A., 1974. Monitoring the vernal advancement and
Landsat-7 Enhanced thematic Mapper Plus (ETM+) and Landsat-8 operational land retrogradation of natural vegetation. NASA/GDFC Type III Final Report, Greenbelt.
imager (OLI) sensors. Remote Sens. 6, 310–329. Royer, A., Charbonneau, L., Brochu, R., Murphy, J.M., Tiellet, P.M., 1987. Radiometric
Lieth, H., 1975. Modeling the primary productivity of the world. In: Lieth, H., Whittaker, comparison of the Landsat-5 TM and MSS sensors. Int. J. Remote Sens. 8, 579–591.
R.H. (Eds.), Primary Productivity in the Biosphere. Springer-Verlag, New York, Running, S.W., Gower, S.T., 1991. Forest-BGC, a general-model of forest ecosystem
237–263. processes for regional applications0.2. Dynamic carbon allocation and nitrogen
Lu, H., Raupach, M.R., McVicar, T.R., Barrett, D.J., 2003. Decomposition of vegetation budgets. Tree Physiol. 9, 147–160.
cover into woody and herbaceous components using AVHRR NDVI time series. Ruprecht, J.K., Schofield, N.J., Crombie, D.S., Stoneman, G.L., 1991. Early hydrological
Remote Sens. Environ. 86, 1–18. response to intense forest thinning in Southwestern Australia. J. Hydrol. 127,
Macfarlane, C., Landman, P., 2002. Final report on establishment of permanent sites for 261–277.

13
C. Macfarlane et al. Remote Sensing Applications: Society and Environment 6 (2017) 1–14

Ruprecht, J.K., Stoneman, G.L., 1993. Water yield issues in the jarrah forest of south- index in Terai forest plantations using fine- and moderate-resolution satellite data.
western Australia. J. Hydrol. 150, 369–391. Int. J. Remote Sens. 35, 7749–7762.
Ryu, Y., Baldocchi, D.D., Kobayashi, H., van Ingen, C., Li, J., Black, T.A., Beringer, J., Turner, D.P., Cohen, W.B., Kennedy, R.E., Fassnacht, K.S., Briggs, J.M., 1999.
Van Gorsel, E., Knohl, A., Law, B.E., Roupsard, O., 2011. Integration of MODIS land Relationships between leaf area index and Landsat TM spectral vegetation indices
and atmosphere products with a coupled-process model to estimate gross primary across three temperate zone sites. Remote Sens. Environ. 70, 52–68.
productivity and evapotranspiration from 1 km to global scales. Glob. Biogeochem. US Geological Survey, 2012. Landsat 7 Science Data Users Handbook. Available online at
Cycles 25, GB4017. 〈http://doi.org/landsathandbook.gsfc.nasa.gov/pdfs/Landsat7_Handbook.pdf〉
Schmidt, M., King, E.A., McVicar, T.R., 2008. A method for operational calibration of US Geological Survey, 2015. Landsat 8 (L8) Data Users Handbook. Available online at
AVHRR reflective time series data. Remote Sens. Environ. 112, 1117–1129. 〈http://landsat.usgs.gov/l8handbook.php〉
Schofield, N.J., Stoneman, G.L., Loh, I.C., 1989. Hydrology of the jarrah forest. In: Dell, Wallace, J.F.,1992. Estimation of leaf area index using remotely sensed image data. A
B. (Ed.), The Jarrah Forest. Kluwer Academic Publishers, Dordrecht, 179–201. report for Alcoa Australia. CSIRO Australia, November 1992.
Schott, J.R., Salvaggio, C., Volchock, W.J., 1988. Radiometric scene normalization using van Leeuwen, W.J.D., Barron, J.O., Marsh, S.E., Herrman, S.M., 2006. Multi-sensor
pseudoinvariant features. Remote Sens. Environ. 26, 1–16. NDVI data continuity: uncertainties and implications for vegetation monitoring
Sellers, P.J., Dickinson, R.E., Randall, D.A., Betts, A.K., Hall, F.G., Berry, J.A., Collatz, applications. Remote Sens. Environ. 100, 67–81.
G.J., Denning, A.S., Mooney, H.A., Nobre, C.A., Sato, N., Field, C.B., Henderson- Wallace, J., 1996. Relationships between satellite image data and leaf area index in jarrah
Sellers, A., 1997. Modeling the exchanges of energy, water, and carbon between forest. A Report to Alcoa Australia. CSIRO Australia.
continents and the atmosphere. Science 275, 502–509. Wang, Y.P., Kowalczyk, E., Leuning, R., Abramowitz, G., Raupach, M.R., Pak, B., van
Silberstein, R.P., Aryal, S.K., Durrant, J., Pearcey, J.M., Braccia, M., Charles, S.P., Gorsel, E., Luhar, A., 2011. Diagnosing errors in a land surface model (CABLE) in
Boniecka, L., Hodgson, G.A., Bari, M.A., Viney, N.R., McFarlane, D.J., 2012. Climate the time and frequency domains. J. Geophys. Res. 116, G01034.
change and runoff in southwestern Australia. J. Hydrol. 475, 441–455. Waring, R.H., 1983. Estimating forest growth and efficiency in relation to canopy leaf
Smettem, K.R.J., Waring, R.H., Callow, J.N., Wilson, M., Mu, Q., 2013. Satellite-derived area. Adv. Ecol. Res. 13, 327–354.
estimates of forest leaf area index in southwest Western Australia are not tightly Watson, D.J., 1958. The dependence of net assimilation rate upon leaf area index. Ann.
coupled to interannual variations in rainfall: implications for groundwater decline in Bot., Lond., New Ser. 22, 37–54.
a drying climate. Glob. Change Biol. 19, 2401–2412. Whitford, K., Colquhoun, I.J., Lang, A.R., Harper, B.M., 1995. Measuring leaf area index
Sokal, R.R., Rohlf, F.J., 1995. Biometry. W.H. Freeman, New York, NY, USA. in a sparse eucalpyt forest: a comparison of estimates from direct measurement,
Song, C., Woodcock, C.E., Seto, K.C., Lenney, M.P., Macomber, S.A., 2001. Classification hemispherical photography, sunlight transmittance and allometric regression. Agric.
and change detection using Landsat TM data: when and how to correct atmospheric For. Meteorol. 74, 237–249.
effects? Remote Sens. Environ. 75, 230–244. Woodcock, C.E., Allen, R., Anderson, M., Belward, A., Bindschadler, R., Cohen, W., Gao,
Soudani, K., François, C., Le Maire, G., Le Dantec, V., Dufrene, E., 2006. Comparative F., Goward, S.N., Helder, D., Helmer, E., Nemani, R., Oreapoulos, L., Schott, J.,
analysis of IKONOS, SPOT, and ETM+ data for leaf area index estimation in Thenkabail, P.S., Vermote, E.F., Vogelmann, J., Wulder, M.A., Wynne, R.H., 2008.
temperate coniferous and deciduous forest stands. Remote Sens. Environ. 102, Free access to Landsat imagery. Science 320, 1011.
161–175. Wulder, M.A., Masek, J.G., Cohen, W.B., Loveland, T.R., Woodcock, C.E., 2012. Opening
Steven, M.D., Malthus, T.J., Baret, F., Xu, H., Chopping, M.J., 2003. Intercalibration of the archive: how free data has enabled the science and monitoring promise of
vegetation indices from different sensor systems. Remote Sens. Environ. 88, Landsat. Remote Sens. Environ. 122, 2–10.
412–422. Xu, D., Guo, X., 2014. Compare NDVI extracted from Landsat 8 imagery with that from
Stoneman, G.L., Schofield, N.J., 1989. Silviculture for water production in jarrah forest of Landsat 7 imagery. Am. J. Remote Sens. 2, 10–14.
Western Australia: an evaluation. For. Ecol. Manag. 27, 273–293. Yoshioka, H., 2003. An isoline-based translation technique of spectral vegetation index
Tiellet, P.M., Barker, J.L., Markham, B.L., Irish, R.R., Fedosejevs, G., Storey, J.C., 2001. using EO-1 hyperion data. IEEE Trans. Geosci. Remote Sens. 41, 1363–1372.
Radiometric cross-calibration of the Landsat-7 ETM+ and Landsat-5 TM sensors Zarate-Valdez, J.L., Whiting, M.L., Lampinen, B.D., Metcalf, S., Ustin, S.L., Brown, P.H.,
based on tandem data sets. Remote Sens. Environ. 78, 39–54. 2012. Prediction of leaf area index in almonds by vegetation indexes. Comput.
Tripathi, P., Patel, N.R., Kushwaha, S.P.S., Dadhwal, V.K., 2014. Upscaling of leaf area Electron. Agric. 85, 24–32.

14