remote sensing
Article

Distinguishing Land Change from Natural Variability
and Uncertainty in Central Mexico with MODIS EVI,
TRMM Precipitation, and MODIS LST Data
Zachary Christman 1, *, John Rogan 2 , J. Ronald Eastman 2,3 and B. L. Turner II 4
1
2
3
4

*

Department of Geography and Environment, Rowan University, Glassboro, NJ 08028, USA
Graduate School of Geography, Clark University, 950 Main Street, Worcester, MA 01610, USA;
jrogan@clarku.edu (J.R.); reastman@clarku.edu (J.R.E.)
Clark Labs, Clark University, 950 Main Street, Worcester, MA 01610, USA
School of Geographical Sciences and Urban Planning, Arizona State University, COOR 5628, Tempe,
AZ 85287-0104, USA; Billie.L.Turner@asu.edu
Correspondence: christmanz@rowan.edu; Tel.: +1-856-256-4810

Academic Editors: Yudong Tian, Ken Harrison, Alfredo R. Huete and Prasad S. Thenkabail
Received: 30 March 2016; Accepted: 2 June 2016; Published: 7 June 2016

Abstract: Precipitation and temperature enact variable influences on vegetation, impacting the
type and condition of land cover, as well as the assessment of change over broad landscapes.
Separating the influence of vegetative variability independent and discrete land cover change remains
a major challenge to landscape change assessments. The heterogeneous Lerma-Chapala-Santiago
watershed of central Mexico exemplifies both natural and anthropogenic forces enacting variability
and change on the landscape. This study employed a time series of Enhanced Vegetation Index (EVI)
composites from the Moderate Resolution Imaging Spectoradiometer (MODIS) for 2001–2007 and
per-pixel multiple linear regressions in order to model changes in EVI as a function of precipitation,
temperature, and elevation. Over the seven-year period, 59.1% of the variability in EVI was explained
by variability in the independent variables, with highest model performance among changing and
heterogeneous land cover types, while intact forest cover demonstrated the greatest resistance to
changes in temperature and precipitation. Model results were compared to an independent change
uncertainty assessment, and selected regional samples of change confusion and natural variability
give insight to common problems afflicting land change analyses.
Keywords: vegetation; variability; Land Use and Land Cover Change; precipitation; temperature;
MODIS; TRMM; EVI; LST

1. Introduction
While the influence on vegetative vigor of climatic variables such as precipitation and temperature
relates to the type and condition of land cover, the isolation of this variability independent of discrete
land cover change is a major challenge of the geographic information science, remote sensing, and land
change science communities. In the heterogeneous and rapidly changing Lerma-Chapala-Santiago
watershed of Mexico [1–4], natural and anthropogenic factors, including droughts, the El Niño
Southern Oscillation, agricultural expansion, and forest loss enact changes in the remotely sensed
detection of vegetative vigor. This study assessed variability using a time series of Enhanced
Vegetation Index (EVI) composites from the Moderate Resolution Imaging Spectoradiometer (MODIS)
for 2001–2007 and per-pixel multiple linear regressions in order to model changes in EVI as a function
of precipitation, temperature, and elevation.
The natural variability of vegetation vigor due to the individual and compound effects of
temperature and precipitation has been investigated in numerous geographical contexts, primarily
Remote Sens. 2016, 8, 478; doi:10.3390/rs8060478

www.mdpi.com/journal/remotesensing

Remote Sens. 2016, 8, 478

2 of 18

through the use of ground-based observation stations or data from the Advanced Very High Resolution
Radiometer (AVHRR) [5]. The need for both inter- and intra-annual analyses for both short-term
assessments and long term monitoring and prediction has been affirmed by many researchers, but
systematic data errors and coarse global datasets have hindered applicability of studies to local and
regional contexts [6–8]. Still, the direct and immediate relationship of increased precipitation to
positive changes in vegetative vigor has been established through the use of remotely sensed data in
mid-latitude semi-arid environments, such as those in central Mexico [9,10]. Because natural variability
may diminish or exaggerate the ability of researchers to monitor discrete changes in land cover with
satellite-derived image products [11], it is important to consider the climatic conditions concurrent
with land change analyses [12] to avoid both errors of omission and commission in land cover change
detection and quantification.
1.1. Vegetation Indices and Vegetative Vigor
Improving the ability to measure the quantity and condition of vegetation from remote sensing
instruments has been an intensive endeavor for more than 40 years. Early efforts by ecological pioneers
such as Jordan [13] and Sellers [14] to relate the quality of light reaching the forest floor to metrics
like leaf area index, have been refined into vegetative indices (e.g., the work of Tucker [15]), which
attempts to maximize the difference between reflected variables to generate metrics of vegetative vigor.
The investigation of phenological variability, including changes to the timing and duration of growing
seasons, has been facilitated through the increased availability of continuous temporal datasets, such
as AVHRR [16,17] and MODIS [18]. Previous studies have indicated changes in the length and onset of
growing seasons around the world, which have been ascribed to interannual variability or long-term
climatic trends [17,19]. These changing intra-annual patterns are ascribable to both climactic variability
and changing land uses [18,19], and regions experiencing these changes are highly susceptible to
misidentifications of land change [20].
The Normalized Difference Vegetation Index (NDVI) has been the most widely used vegetation
index, and its broad utility across different land covers and remote sensing instruments has led to
tremendous gains in the understanding of regional and global land cover [21,22]. Over tropical and
subtropical regions, the NDVI ratio, using satellite image bands in the red and near-infrared ranges,
is subject to saturation in cases of high vegetative vigor and sensitivity to atmospheric conditions [22].
Numerous efforts to quantify the relationship between reflected radiation and ecological variables
such as leaf-area index and vegetation cover [5,14] have led to improvements upon the NDVI [22]. The
Enhanced Vegetation Index (EVI) (designed to minimize haze effects on red reflectance and expand
radiometric sensitivity to varying amounts of vegetation, by incorporating blue band and correction
constant values) has been systematically evaluated and compared to NDVI in tropical and subtropical
environments, and EVI has proven robust to saturation in regions of dense vegetation, while remaining
sensitive to variations in canopy cover [23,24]. Several studies have compared the effectiveness of
NDVI and EVI in tropical [25] and semi-arid environments [26], noting the close relationship between
the two indices. In a study of the effect of sub-pixel mixing of land cover on the values of vegetation
indices, Liu and Kafatos [27] found that EVI varied less than NDVI in areas of mixed land cover,
making it especially suitable in heterogeneous landscapes. Therefore, the MODIS monthly 1-kilometer
EVI composite product was chosen to represent the vegetative vigor of land cover in relation to climate
variables of temperature and precipitation.
1.2. Linear Modeling in Time Series Analysis
Linear modeling, in which predictor equations are calculated based on the correlations of ordinary
least squares regression, has proven useful in identifying and quantifying the relationships among
time series of bioclimatic variables [28]. The use of linear modeling in the context of time series
analysis identifies the relationship among temporal sequences of imagery, rather than individual

Remote Sens. 2016, 8, 478

3 of 18

images. Rather than a single slope and intercept to define the terms of the relationship, the relationship
of the dependent variable to the independent variables is calculated on a per-pixel basis.
The relationship of temperature and precipitation with variation in vegetation indices has been
evaluated using AVHRR [6,29] and ground-based data [12]. In almost every global ecosystem, there is a
positive relationship between increased precipitation and increased vegetation index values, especially
in water-stressed regions [6,30,31]. Generally, in both semi-arid and sub-tropical regions, a negative
relationship between temperature and EVI indicates vegetation vigor is limited by water availability
and especially soil moisture [7,8,32]. In the investigation of time series relationships, a change in an
independent variable may occur concurrently with a change in the dependent variable. Alternately,
a lagged effect implies that a change in an independent variable may precede a change in the dependent
variable or that changes in the dependent variable were ascribable to the compound effects of the
preceding changes in the independent variables. Wang and colleagues have demonstrated that while
there are lagged effects of precipitation and temperature upon vegetation variation, the highest
correlation occurred within one or two biweekly periods, or approximately one month [31]. Further,
Braswell and colleagues note that there are both significant instantaneous and lagged relationships
among vegetation and temperature and that these vary across ecosystems [7]. While lagged links
between EVI and climate variables may represent subtle ongoing processes of change [33], they are not
effective for the identification of discrete categorical change due to the longer timelines in which these
changes occur. As the impacts of precipitation on the natural and rain fed agricultural landscapes of
this region would have a quicker impact than the temporal resolution of the data, this study addresses
the instantaneous changes in temperature and precipitation that occur concurrently with changes
in EVI in order to relate the changes in vegetative vigor to the identification of discrete change in
categorical land cover assessments.
1.3. Study Area
This research was conducted in the Lerma-Chapala-Santiago watershed (Figure 1) of central
Mexico (19–23.5˝ N, 99–105.5˝ W), on Mexico’s Central Mesa, bordered to the south by the Transverse
Neovolcanic Mountain Range [34]. Encompassing nearly 136,000 square kilometers, the watershed
spans more than 4000 meters in elevation, including the drainages of the Lerma River (Río Lerma) and
the Santiago River (Río Grande de Santiago) [35].
This region includes a wide variety of climatic zones and natural vegetation types, due to its
diverse topography and geomorphology [34]: semi-arid, dry-steppe, humid sub-tropical, pine-oak
forests, and others. As a region situated within the tropics, cyclical patterns of natural variability, such
as the El Niño Southern Oscillation, have a direct and immediate effect on precipitation regimes [36].
The phase of El Niño brings cooler temperatures and increased precipitation, while the complementary
phase of La Niña is associated with dryer, hotter conditions [36]. During this study, a minor La Niña
pattern was identified in 2000–2001 and 2007, and an El Niño pattern has been identified in the winters
of 2002–2003, 2004–2005, and 2006–2007. Additionally, the region is subject to the potential influence of
both Atlantic (Caribbean) and Pacific hurricanes, though its elevation and relative topography limit
the influence of these forces largely to changes in precipitation [36,37].
Patterns indicating widespread drying phenomena have been noted across southwestern North
America [38], and there are indications of a major drought in Mexico in the early twenty-first
century [39]. Across central Mexico, there are numerous unprotected water bodies that qualify
as threatened wetland ecosystems under international conventions [40], and large lakes in this study
area, including Chapala and Cuitzeo, have seen volumes decline at an unprecedented rate over recent
decades [35,41]. Drought events have a major impact on the Mexican agricultural sector and the
communities who depend on it [42,43], and the linear patterns of increased drying and decreased
precipitation are shaping a variety of industries [44].

Remote Sens. 2016, 8, 478
Remote Sens. 2016, 8, 478

4 of 18

4 of 17

Figure
1. Location
and
ElevationofofLerma-Chapala-Santiago
Lerma-Chapala-Santiago Watershed
andand
fourfour
selected
sample
Figure
1. Location
and
Elevation
Watershed
selected
sample
regions (a–d), with location in country of Mexico.
regions (a–d), with location in country of Mexico.

The natural heterogeneity of the region is exemplified by a wide variety of plant species. There
Themore
natural
the region
is exemplified
by ainwide
variety of
of Mexico
plant species.
are
thanheterogeneity
40 species and of
recognized
variations
of pine trees
the Republic
[45], andThere
the Transvolcanic
range
Central Mesa
are recognized
as centers
of pine species
diversity
in the
are more
than 40 species
andand
recognized
variations
of pine trees
in the Republic
of Mexico
[45], and
Mexico
[46].
The
Monarch
Butterfly
Biosphere
Reserve,
on
the
edge
of
the
watershed
near
the
borders
Transvolcanic range and Central Mesa are recognized as centers of pine species diversity in Mexico [46].
of the states
of Michoacán
and Mexico,
represents
a unique
realized
niche critical
to borders
the annual
The Monarch
Butterfly
Biosphere
Reserve,
on the edge
of the
watershed
near the
ofcycle
the states
of the global population of the species [47].
of Michoacán and Mexico, represents a unique realized niche critical to the annual cycle of the global
The broad range of natural variation in land cover types and the apparent fluctuation in vegetation
population of the species [47].
vigor across the region is exacerbated by human activity. In addition to extensive heterogeneity in its
The broad range of natural variation in land cover types and the apparent fluctuation in vegetation
natural cover, the Lerma-Chapala-Santiago watershed includes several major urban centers, including
vigorthe
across
region is exacerbated
by humanToluca
activity.
In addition
to extensive
cities the
of Guadalajara
(1.7 million inhabitants),
(~750,000
inhabitants),
Queretaro heterogeneity
(~730,000
in itsinhabitants),
natural cover,
the
Lerma-Chapala-Santiago
watershed
includes
several
major
urban
centers,
Aguascalientes (~700,000 inhabitants), Morelia (~640,000 inhabitants), and Guanajuato
including
the
cities
of
Guadalajara
(1.7
million
inhabitants),
Toluca
(~750,000
inhabitants),
Queretaro
(~480,000 inhabitants) (INEGI, 2005). With many regional urban centers within the watershed and
Mexicoinhabitants),
City just beyond
its limits, the expansion
of inhabitants),
the built environment
exerts
a stronginhabitants),
and constant and
(~730,000
Aguascalientes
(~700,000
Morelia
(~640,000
pressure (~480,000
on natural inhabitants)
land covers [48,49]
and 2005).
water resources
[35,50,51].
Processes
urbanization,
Guanajuato
(INEGI,
With many
regional
urbanofcenters
within the
extensive
and
intensive
agriculture,
timber
harvesting,
and
mining
fragment
the
landscape,
creating
watershed and Mexico City just beyond its limits, the expansion of the built environment exerts a
a motley patchwork of human and natural land covers across most of the watershed [3,52,53].
strong
and constant pressure on natural land covers [48,49] and water resources [35,50,51]. Processes
Additionally, this agriculturally-intensive region has seen a trend toward the consolidation of
of urbanization, extensive and intensive agriculture, timber harvesting, and mining fragment the
smallholders’ plots into more commercial and mechanized agricultural systems [54,55]. The
landscape,
creating a motley patchwork of human and natural land covers across most of the
watershed is a major producer of corn for the region and the country, and increased investment in
watershed
[3,52,53].
agriculturally-intensive
region hasinseen
a trendstate
toward
tequila productionAdditionally,
has led to thethis
conversion
of fields to agave cultivation
the western
of the
consolidation
of
smallholders’
plots
into
more
commercial
and
mechanized
agricultural
systems
[54,55].
Jalisco [56,57].

The watershed is a major producer of corn for the region and the country, and increased investment
in tequila production has led to the conversion of fields to agave cultivation in the western state of
Jalisco [56,57].

Remote Sens. 2016, 8, 478
Remote Sens. 2016, 8, 478

5 of 18

5 of 17

Many landscapes in the region comprise a heterogeneous
heterogeneous mixture
mixture of
of multiple
multiple covers.
covers. In the
classification
scheme
of
the
International
Geosphere-Biosphere
Program
Data
and
Information
System
classification scheme of the International Geosphere-Biosphere Program Data and Information
(IGBP-DIS)
DISCover
land cover
product,
a thematic
class of “Cropland/Natural
Vegetation
Mosaic”
System (IGBP-DIS)
DISCover
land
cover product,
a thematic
class of “Cropland/Natural
Vegetation
was
included
to accommodate
cover types
that
did that
not represent
pure endmembers
of eitherofnatural
Mosaic”
was included
to accommodate
cover
types
did not represent
pure endmembers
either
or
anthropogenic
land
covers
[58,59].
This
“Mosaic”
cover
is
especially
appropriate
in
the
diverse
natural or anthropogenic land covers [58,59]. This “Mosaic” cover is especially appropriate in the
Lerma-Chapala-Santiago
watershed
as a buffer
between
natural
land
covers,
as
diverse
Lerma-Chapala-Santiago
watershed
as class
a buffer
class cropland
between and
cropland
and
natural
land
well
as serving
as serving
an indicator
forest fragmentation
or harvesting
[58,60].
covers,
as well as
as anof
indicator
of forest fragmentation
or harvesting
[58,60].
2. Materials
Materialsand
andMethods
Methods
2.
Data for
for this
index and
products from
Data
this project
project included
included vegetation
vegetation index
and land
land surface
surface temperature
temperature products
from
MODIS,
average
daily
precipitation
from
TRMM,
a
DEM
from
SRTM,
and
maps
of land
MODIS, average daily precipitation from TRMM, a DEM from SRTM, and maps of
land cover
cover
persistence
Methods explicated
explicated below
below included
included the
the extraction
extraction of
of summary
persistence and
and change
change [1].
[1]. Methods
summary statistics
statistics
from
the
time
series
of
dependent
and
independent
variables
and
a
multiple
linear
regression,
from the time series of dependent and independent variables and a multiple linear regression,
performed
seven-year
series
and and
each each
annualannual
subseries
independently
(summarized
patterns
performed over
overthe
the
seven-year
series
subseries
independently
(summarized
of EVI, precipitation,
and temperature
shown in
Figure
Results
of the of
regression
modelmodel
were
patterns
of EVI, precipitation,
and temperature
shown
in 2).
Figure
2). Results
the regression
compared
to
a
previous
land
cover
and
uncertainty
analysis
to
explain
the
effect
of
variability
upon
were compared to a previous land cover and uncertainty analysis to explain the effect of variability
land change
error and
[1]. [1].
upon
land change
errorconfusion
and confusion

Figure
Monthly
Enhanced
Vegetation
Index (EVI),
and Temperature,
2001–2007,
Figure 2.2.Average
Average
Monthly
Enhanced
Vegetation
IndexPrecipitation,
(EVI), Precipitation,
and Temperature,
for
entire the
watershed.watershed.
2001–2007,
forLerma-Chapala-Santiago
entire the Lerma-Chapala-Santiago

2.1.
2.1. Indices
Indices of
of Vegetative
Vegetative Vigor
Vigor
Enhanced
Index (EVI)
(EVI) image
image composites
composites produced
produced by
by MODIS
MODIS product
product MOD13A3
MOD13A3
Enhanced Vegetation
Vegetation Index
were
used
for
this
study.
This
product
is
a
monthly
maximum
value
composite
of
observations
were used for this study. This product is a monthly maximum value composite of observations from
from
the
MODISinstrument
instrument
Terra
platform,
a spatial
resolution
of 1pixel
km and
percontinuous
pixel and
the MODIS
onon
thethe
Terra
platform,
with awith
spatial
resolution
of 1 km per
continuous
coverage
over of
thestudy,
period
of study, 2001–2007.
As this
imagery
has undergone
geometric
coverage over
the period
2001–2007.
As this imagery
has
undergone
geometric correction
correction
and
compositing
by
the
producers
to
accommodate
periods
of
impaired
observations
[61],
and compositing by the producers to accommodate periods of impaired observations [61], no further
no
further was
correction
was performed
in this study.
correction
performed
in this study.
2.2. Temperature, Precipitation and Elevation Data
Land surface temperature data were derived from the MODIS Terra product MOD11C3, which
are distributed as monthly composites of 0.05 degrees per pixel. Data units were transformed from

Remote Sens. 2016, 8, 478

6 of 18

2.2. Temperature, Precipitation and Elevation Data
Land surface temperature data were derived from the MODIS Terra product MOD11C3, which
are distributed as monthly composites of 0.05 degrees per pixel. Data units were transformed from
Remote Sens. 2016, 8, 478
of 17
Kelvin
to degrees Celsius to maximize interpretability of the results, as negative temperature6 values
were not present in the sample.
Kelvin to degrees Celsius to maximize interpretability of the results, as negative temperature values
Precipitation data were derived from the Tropical Rainfall Measuring Mission product 3B43,
were not present in the sample.
which are distributed as monthly composites of average hourly precipitation at 0.25˝ per pixel [62].
Precipitation data were derived from the Tropical Rainfall Measuring Mission product 3B43,
Precipitation
data values
rescaled
from mm/hour
mm/day
to maximize
interpretability
which are distributed
aswere
monthly
composites
of average to
hourly
precipitation
at 0.25°
per pixel [62]. of
thePrecipitation
results.
data values were rescaled from mm/hour to mm/day to maximize interpretability of
To
account
for the effects of elevation upon temperature and precipitation [63,64], a digital
the results.
elevation
model
(DEM)
from
the of
Shuttle
Radar
Topography
Mission
was included[63,64],
in the aregression
To account for the
effects
elevation
upon
temperature
and precipitation
digital
analysis
as amodel
constant.
series
identical
SRTM
images Mission
was used
inincluded
the regression
to mimic a
elevation
(DEM)Afrom
theofShuttle
Radar
Topography
was
in the regression
temporal
of theAstatic
A random
factor of
less
than
meter wasto
added
analysissequence
as a constant.
seriesvariable.
of identical
SRTM images
was
used
in one
the regression
mimictoaall
SRTM
images
to avoidofthe
of a singular
correlation
matrix
this series.
temporal
sequence
thecalculation
static variable.
A random
factor of less
thanamong
one meter
was added to all
SRTM
images geometric
to avoid thecompatibility
calculation of awith
singular
among
this series.
To ensure
the correlation
series of matrix
vegetation
indices,
all temperature,
To ensure
compatibility
with the series
of vegetation
indices, all temperature,
precipitation,
and geometric
elevation images
were geometrically
reprojected
and downsampled
to 1 km per
precipitation,
and
elevation
images
were geometrically
and downsampled
to 1 km per
pixel
resolution in
the
sinusoidal
references
system usingreprojected
a bilinear calculation
[65].
pixel resolution in the sinusoidal references system using a bilinear calculation [65].
2.3. Land Cover Dynamics
2.3. Land Cover Dynamics
Maps of eight natural and anthropogenic land cover categories for 2001 and 2007 were
Maps ofto eight
natural
and
anthropogenic
land cover
categories
forthe
2001
and 2007
were[1]
crosstabulated
identify
regions
of persistent
and changing
land
cover over
seven-year
period
crosstabulated
to
identify
regions
of
persistent
and
changing
land
cover
over
the
seven-year
period
[1]
(Figure 3). Land cover classes were based on a systematic roadside survey of the Lerma-Chapala(Figure
3).
Land
cover
classes
were
based
on
a
systematic
roadside
survey
of
the
Lerma-ChapalaSantiago during January and July 2007 and included the following thematic categories: forest; shrub;
Santiago
during
Januarybuilt
and (including
July 2007 and
included
the following
thematic
categories:
forest;
shrub;
grass;
sparse
vegetation;
urban
and rural
settlements);
crop
(including
mechanized
grass; sparse vegetation; built (including urban and rural settlements); crop (including mechanized
and smallholder agriculture); and mosaic (defined by a sub-pixel mixture of cropland and natural
and smallholder agriculture); and mosaic (defined by a sub-pixel mixture of cropland and natural
vegetation). Land cover maps were generated for 2001 and 2007 using a combined pseudo-invariant
vegetation). Land cover maps were generated for 2001 and 2007 using a combined pseudo-invariant
calibration database and the Mahalanobis distance classifier in the Idrisi Taiga software package [1].
calibration database and the Mahalanobis distance classifier in the Idrisi Taiga software package [1].
Overall map accuracy for the both maps was approximately 90%.
Overall map accuracy for the both maps was approximately 90%.

(a)

(b)

Figure 3. Selected regions of high confidence land cover persistence (a) and change from 2001 to 2007
Figure 3. Selected regions of high confidence land cover persistence (a) and change from 2001 to 2007
(b) across the Lerma-Chapala-Santiago watershed.
(b) across the Lerma-Chapala-Santiago watershed.

Regions of changing and persistent land cover were further refined using the Confusion Index
Regions
of changing
and persistent
land cover
were further
Index
metric
derived
from soft-classified
Mahalanobis
typicality
imagesrefined
for bothusing
dates the
[1]. Confusion
The Confusion
metric
derived
from soft-classified
typicality
for both
dates
[1].
Confusion
Index
(CI) quantifies
the potentialMahalanobis
for error in land
change images
assessments
based
upon
theThe
uncertainty
inherent to each map used for the change assessment, identifying potential cases in which
classification error in one or both maps may contribute to the misidentification of landscape change.
CI values range from 0 to 1, in which low values represent a low likelihood of spurious change
assessments and higher values indicate apparent change (or persistence) in the imagery may not

Remote Sens. 2016, 8, 478

7 of 18

Index (CI) quantifies the potential for error in land change assessments based upon the uncertainty
inherent to each map used for the change assessment, identifying potential cases in which classification
error in
one
or8,both
Remote
Sens.
2016,
478 maps may contribute to the misidentification of landscape change. CI values
7 of 17
range from 0 to 1, in which low values represent a low likelihood of spurious change assessments and
match
landscape
For this(or
study,
land cover
persistence
and change
assessments
were
higher true
values
indicatedynamics.
apparent change
persistence)
in the
imagery may
not match
true landscape
thresholded
at athis
level
below
0.3,cover
based
on the results
a sensitivity
analysis,
to thresholded
isolate the transitions
dynamics. For
study,
land
persistence
andof
change
assessments
were
at a level
in
which
was
confidence.
below
0.3,there
based
on the
the highest
results of
a sensitivity analysis, to isolate the transitions in which there was the
highest confidence.
2.4. Multiple Linear Regression
2.4. Multiple Linear Regression
Multiple linear regression was used to model the series of monthly EVI composites as a
Multiple
linear against
regression
used to model
the series
of monthly and
EVI composites
asThrough
a dependent
dependent
variable
thewas
independent
variables
of precipitation
temperature.
the
variable
against
independent
of precipitation
and temperature.
Through
theslope
use ofand
the
use
of the
Earththe
Trends
Modelervariables
of the Idrisi
Taiga software
package [66],
per-pixel
Earth Trends
Modeler
of the Idrisi
package [66],coefficient
per-pixel slope
and interceptfor
terms
intercept
terms
were calculated
asTaiga
well software
as a spatially-explicit
of determination
the
were calculated
as wellModels
as a spatially-explicit
coefficient
determination
modeledasrelationship.
modeled
relationship.
were generated
using theofseven-year
seriesfor
ofthe
84 months
well seven
Models were annual
generated
using
themonths/year.
seven-year series
84 months
well sevenseries
independent
independent
series
of 12
Modelofresults
of theasseven-year
(R2) are annual
shown
series
of 12
in
Figure
4. months/year. Model results of the seven-year series (R2 ) are shown in Figure 4.

2
Figure
4. Coefficient
Figure 4.
Coefficient ofof determination
determination (R
(R)2 ) ofof modeled
modeled variability
variability in
in EVI
EVI with
with temperature,
temperature,
precipitation,
precipitation, and
and elevation
elevation as
as independent
independent variables
variables in
in 2001–2007
2001–2007 series
series across
across the
the Lerma-ChapalaLerma-ChapalaSantiago
Santiago watershed.
watershed.

2.5.
2.5. Statistics
Statistics of
of Variability
Variability and
and Model
Model Performance
Performance
Summary
Summary statistics
statistics of
of the
the average
average monthly
monthly of
of EVI,
EVI, temperature,
temperature, and
and precipitation
precipitation and
and the
the annual
annual
variability
of
each
of
these
variables
were
calculated
over
the
entire
study
area.
The
mean
and
variability of each of these variables were calculated over the entire study area. The mean and
standard
standard
of the monthly
spatial averages
was calculated
to indicate
the temporal
variability
deviationdeviation
of the monthly
spatial averages
was calculated
to indicate
the temporal
variability
of each
of
each
variable
for
each
location
through
time
on
an
annual
basis.
Additionally,
the
spatial
variable for each location through time on an annual basis. Additionally, the spatial population
population
standardwas
deviation
was from
calculated
from the
annual
mean EVI
images tothe
illustrate
the
standard deviation
calculated
the annual
mean
EVI images
to illustrate
variability
variability
across
space.
across space.
The spatially-explicit model equation terms and the coefficient of determination were averaged
over the scale of the entire watershed as well as the isolated regions of persistent and changing land
cover, for the seven-year series and the annual series relationships.
2.6. Selected Sample Regions
Four selected sample regions were isolated based on model performance and Confusion Index
(CI) values to illustrate potential cases in which land cover change and persistence may or may not

Remote Sens. 2016, 8, 478

8 of 18

The spatially-explicit model equation terms and the coefficient of determination were averaged
over the scale of the entire watershed as well as the isolated regions of persistent and changing land
cover, for the seven-year series and the annual series relationships.
2.6. Selected Sample Regions
Four selected sample regions were isolated based on model performance and Confusion Index
Remote
Sens. 2016,
8, 478
8 of be
17
(CI)
values
to illustrate
potential cases in which land cover change and persistence may or may not
conflated with natural variability. Regions were selected quantitatively through the logical intersection
2 values,
intersection
veryand
lowvery
(<0.4)high
and (>0.7)
very high
(>0.7) CI,
whichfrom
ranged
0 tomodel
1, andRmodel
R2 values,
of very low of
(<0.4)
CI, which
ranged
0 tofrom
1, and
which
which
0 to
0.84.with
Areas
with
CI comprised
of the landscape,
while
areas
with
rangedranged
from 0 from
to 0.84.
Areas
low
CI low
comprised
22.0% of22.0%
the landscape,
while areas
with
high
CI
2
2
high
CI comprised
of the landscape.
Areaslow
with
R values
comprised
the landscape
comprised
49.0% of49.0%
the landscape.
Areas with
R low
values
comprised
7.7% 7.7%
of theoflandscape
and
and
R2 values
comprised
the landscape.
The intersection
yielded
thezones
four from
zoneswhich
from
highhigh
R2 values
comprised
15.3%15.3%
of theoflandscape.
The intersection
yielded
the four
which
samples
were
selected,
with
their
proportional
area
of
the
entire
Lerma-Chapala-Santiago
samples were selected, with their proportional area of the entire Lerma-Chapala-Santiago watershed:
2
2 andlow
watershed:
lowlow
CI, R
2.2%;
R2 CI,
and3.8%;
high high
CI, 3.8%;
high
R2 CI,
and3.6%;
low CI,
3.6%;
highhigh
R2 and
low R2 and low
low RCI,and
2.2%;
high
R2 and
low
high
R2 and
CI,
high
CI,
6.7%.
A
region
of
approximately
21.4
square
kilometers
(corresponding
to
25
1-kilometer
6.7%. A region of approximately 21.4 square kilometers (corresponding to 25 1-kilometer MODIS
MODISorpixels
or 100pixels)
500-mwas
pixels)
was selected
to represent
each sample.
Land
cover composition
pixels
100 500-m
selected
to represent
each sample.
Land cover
composition
in 2001,
2
2
in
2001,
2007,
Confusion
Index,
Model
R
,
Average
EVI,
and
Annual
EVI
profiles
for
2001
and
2007
2007, Confusion Index, Model R , Average EVI, and Annual EVI profiles for 2001 and 2007 were
were
extracted
for
each
sample
(illustrated
in
Figure
5).
extracted for each sample (illustrated in Figure 5).

Figure 5.
Selected samples
samples of
of land
including land
land cover
cover 2001
2001 and
Figure
5. Selected
land change
change scenarios,
scenarios, including
and 2007,
2007, Confusion
Confusion
Index,
Model
Performance,
Average
EVI,
and
Annual
profiles.
(a)
Intact
forest;
(b)
Irrigated
agriculture;
Index, Model Performance, Average EVI, and Annual profiles. (a) Intact forest; (b) Irrigated
(c)
Seasonal (c)
agriculture;
Shrub and
agriculture;
Seasonal (d)
agriculture;
(d)Mosaic.
Shrub and Mosaic.

3. Results
Results
3.
3.1. Average
Average Characteristics
Characteristics of
3.1.
of Dependent
Dependent and
and Independent
Independent Variables
Variables by
by Series
Series
Over the
Over
the seven-year
seven-year period,
period, the
the average
average EVI
EVI value
value across
across the
the entire
entire Lerma-Chapala-Santiago
Lerma-Chapala-Santiago
watershed
was
0.270,
with
a
temporal
standard
deviation
of
0.100
and
a
spatial
standard deviation
deviation of
of
watershed was 0.270, with a temporal standard deviation of 0.100 and a spatial standard
0.408, demonstrating
demonstrating the
the range
range of
Annual average
average EVI
EVI values
0.408,
of variability
variability of
of this
this landscape.
landscape. Annual
values ranged
ranged

from a low of 0.255 in 2001 to a high of 0.284 in 2004. The temporal standard deviation of the sevenyear series was 0.100, which varied on an annual basis from a low of 0.096 in 2006 to a high of 0.105
in 2007. The spatial standard deviation calculated of the averages of the seven-year series was 0.058,
which varied in the annual series from 0.59 (2006) to 0.062 (2004 and 2005).
Average daily precipitation over the entire series was 1.793 mm/day, with a standard deviation

Remote Sens. 2016, 8, 478

9 of 18

from a low of 0.255 in 2001 to a high of 0.284 in 2004. The temporal standard deviation of the seven-year
series was 0.100, which varied on an annual basis from a low of 0.096 in 2006 to a high of 0.105 in 2007.
The spatial standard deviation calculated of the averages of the seven-year series was 0.058, which
varied in the annual series from 0.59 (2006) to 0.062 (2004 and 2005).
Average daily precipitation over the entire series was 1.793 mm/day, with a standard deviation of
2.158 mm/day. On an annual basis, average daily precipitation ranged from a high of 2.579 mm/day in
2004 to a low of 0.769 mm/day in 2005. Average daytime temperature over the seven-year series was
30.357 degrees C, with a standard deviation of 5.980 degrees C. Annually, average daily temperature
ranged from a high of 31.440 degrees C in 2001 to a low of 28.457 degrees C in 2004. In the SRTM
DEM, the elevation of the Lerma-Chapala-Santiago watershed ranges from 4360 m to sea level, with an
average elevation of 1903.9 m. Table 1 details summary statistics of EVI, precipitation, and temperature,
and the profiles of each of these variables is illustrated in Figure 2.
Table 1. Summary statistics of time series of dependent and independent variables across entire
Lerma-Chapala-Santiago watershed.
EVI
Series Year(s)
2001–2007
2001
2002
2003
2004
2005
2006
2007

Mean
0.270
0.255
0.270
0.282
0.284
0.260
0.267
0.276

Precipitation

Standard Deviation
Temporal

Spatial

0.100
0.102
0.101
0.108
0.109
0.098
0.096
0.105

0.058
0.060
0.061
0.060
0.062
0.062
0.059
0.060

Temperature

Mean

Standard
Deviation

Mean

Standard
Deviation

1.793
1.997
2.134
2.433
2.579
0.769
1.227
1.414

2.158
2.291
2.125
3.033
2.810
0.977
1.302
1.644

30.357
31.440
30.309
30.145
28.457
31.386
30.633
30.128

5.980
5.878
6.321
6.519
5.525
6.280
6.713
5.589

3.2. Average Characteristics of Dependent and Independent Series by Land Cover Type and Dynamics
Average EVI values across the sample set of persistent land cover type varied widely. Water, with
little vegetative response, had a mean EVI value of 0.022 across the seven-year series. The highest mean
EVI value was in the forest class, with 0.304, and the lowest mean EVI was in the built class, at 0.151.
Agricultural classes of crop and mosaic had average EVI values of 0.296 and 0.298, respectively.
Precipitation falling on each land cover type varied accordingly with its distribution across the
watershed, ranging from a high of 1.903 mm/day for the mosaic class to a low of 1.395 mm/day for
the grass class. Similarly, average daytime temperature ranged from lows of 22.7 degrees C across the
water cover and 25.5 degrees C across the forest cover to 34.4 across sparsely vegetated land cover.
When averaged over each land cover class, EVI generally exhibited a positive relationship to
elevation (R2 = 0.41), supported by the high average EVI and high elevation of the forest class. Average
temperature and elevation over each land cover showed a weak negative correlation (R2 = 0.03), and
there was no demonstrable relationship between precipitation and elevation over the averaged extent
of each land cover.
For all persistent land covers, the mean EVI for the seven-year series was 0.281, with an average
precipitation of 1.793 mm/day and temperature of 30.3 degrees C. For any land cover that confidently
experienced a categorical shift [1], the mean EVI for the seven-year series was 0.260, with average
precipitation of 1.787 mm/day and temperature of 30.9 degrees C. Table 2 details summary statistics of
average EVI, temperature, and precipitation by land cover type and dynamic over the seven-year series.

Remote Sens. 2016, 8, 478

10 of 18

Table 2. Summary statistics of average EVI, temperature, and precipitation by land cover type and
dynamics, 2001–2007.
Land Cover
Dynamic *

Land
Cover Type

EVI
(mean)

Precipitation
(mean, mm/day)

Temperature
(mean, degrees C)

Elevation
(mean, meters)

Water
Forest
Shrub
Grass
Sparse
Built
Crop
Mosaic

0.022
0.304
0.234
0.204
0.186
0.151
0.296
0.289

1.638
1.799
1.532
1.395
1.587
1.809
1.870
1.903

22.732
25.467
31.975
32.025
34.401
32.578
31.367
30.336

1518.002
2360.186
1918.243
2143.839
1424.184
1694.716
1825.303
1881.013

All Land **

0.281

1.793

30.332

1934.101

All Changes

0.260

1.787

30.903

1975.454

Persistence

Change

* Land Cover Dynamic refers to the state of change or persistence in the comparison of the 2001 and 2007 land
cover maps; ** The “All Land” Land cover type aggregates all land cover types except Water.

3.3. Model Performance by Series
The model indicated that 59.1% of the variability in the dependent variable (EVI) was explained
by precipitation, temperature, and elevation over the seven-year series from 2001 to 2007 (Figure 4).
On an annual basis, EVI exhibited a varying relationship with the independent variables. From
2001 to 2004, model performance across the watershed was high, with 71.5% (2001) to 79.1% (2002)
of the variability in EVI explained by precipitation, temperature, and elevation. In 2005 the model
relationship decreased, with only 52.5% of the variability in EVI explained by independent variables.
In 2006, model performance was highest among all years studied (R2 = 0.798), and in 2007, model
performance decreased slightly, explaining only 68.1% of the variability in EVI. Table 3 summarizes the
annual and series model performances, as well as the slope and intercept coefficients for each model.
Table 3. Average coefficients of determination, slope, and intercept terms for each model across entire
study region.
Series
Year(s)

R2

Slope1
(Precipitation)

Slope2
(Temperature)

Slope3
(Elevation)

Intercept

α

2001–2007
2001
2002
2003
2004
2005
2006
2007

0.5910
0.7148
0.7907
0.7739
0.7246
0.5247
0.7976
0.6809

2.955 ˆ 10´2
2.766 ˆ 10´2
3.837 ˆ 10´2
2.779 ˆ 10´2
3.004 ˆ 10´2
5.564 ˆ 10´2
5.655 ˆ 10´2
4.468 ˆ 10´2

´5.592 ˆ 10´3
´7.549 ˆ 10´3
´3.590 ˆ 10´3
´4.150 ˆ 10´3
´4.079 ˆ 10´3
´5.330 ˆ 10´3
´4.566 ˆ 10´3
´6.139 ˆ 10´3

3.300 ˆ 10´5
1.440 ˆ 10´4
1.740 ˆ 10´4
2.500 ˆ 10´4
´4.500 ˆ 10´5
4.400 ˆ 10´5
´1.180 ˆ 10´4
1.070 ˆ 10´4

0.3376
0.163
0.722
´0.107
0.355
0.270
0.603
0.194

<0.001
<0.001
<0.001
<0.001
<0.001
<0.05
<0.001
<0.01

3.4. Model Performance by Land Cover Type and Dynamics
Across the seven-year series, the average performance of the per-pixel multiple linear regression
varied widely by land cover type. Table 4 summarizes the performance of the seven-year per-pixel
multiple linear regression model by land cover type and dynamics. In a sample of confidently persistent
land covers, on a per-class basis, from 45.9% to 68.0% of the variability in EVI values was explained
by independent variables temperature, precipitation, and elevation. Persistent forested landscapes
showed the lowest modeled relationship, with only 45.9% of the variability in EVI values explained
by the independent variables. The highest model performance was in the mosaic class, explaining
68.0% of the variability in EVI values. Water was excluded from aggregations of land cover, due to
the very low relationship among its vegetative vigor (ascribable mostly to changing water levels and

Remote Sens. 2016, 8, 478

11 of 18

eutrophication), temperature, and precipitation. For any pixel that experienced a categorical change in
land cover type over the seven-year period, 64.5% of the variability in EVI values was explained by
the independent variables. Of the sample of persistent land-based classes, the model explained 56.9%
of the variability in EVI values.
Table 4. Average coefficients of determination, slope, and intercept terms for 2001–2007 series model
by land cover type and dynamics.
Land Cover
Dynamic *

Persistence

Change

Land
Cover Type

R2

Slope1
(Precipitation)

Slope2
(Temperature)

Slope3
(Elevation)

Intercept
Unmodified

α Unmodified

Water
Forest
Shrub
Grass
Sparse
Built
Crop
Mosaic

0.0361
0.4587
0.6053
0.5583
0.5666
0.4423
0.5477
0.6796

´6.010 ˆ 10´4
1.586 ˆ 10´2
3.629 ˆ 10´2
2.808 ˆ 10´2
1.907 ˆ 10´2
1.045 ˆ 10´2
3.194 ˆ 10´2
3.590 ˆ 10´2

´1.558 ˆ 10´3
´2.054 ˆ 10´3
´4.615 ˆ 10´3
´3.437 ˆ 10´3
´1.940 ˆ 10´3
´7.150 ˆ 10´4
´5.415 ˆ 10´3
´7.773 ˆ 10´3

9.030 ˆ 10´4
1.390 ˆ 10´4
´1.060 ˆ 10´4
´3.740 ˆ 10´4
5.250 ˆ 10´4
´1.467 ˆ 10´3
2.010 ˆ 10´4
´1.010 ˆ 10´4

´1.313
7.087 ˆ 10´2
0.574
0.952
´0.559
2.615
´1.763 ˆ 10´2
0.839

<0.4
<0.001
<0.001
<0.001
<0.001
<0.001
<0.001
<0.001

All Land **

0.5690

3.057 ˆ 10´2

´5.097 ˆ 10´3

4.600 ˆ 10´5

0.322

<0.001

0.6453

10´2

10´3

´2.820 ˆ 10´4

1.060

<0.001

All Changes

3.425 ˆ

´6.162 ˆ

* Land Cover Dynamic refers to the state of persistence or change in the comparison of the 2001 and 2007
land cover maps, controlled for uncertainty; ** The “All Land” cover type aggregates all land cover types
except Water.

3.5. Model Performance, Confusion, and Land Change in Selected Samples
Four sample regions were selected to evaluate the cases in which high and low Confusion Index
values intersected with high and low coefficients of determination in the seven-year model (Figure 5).
‚

‚

‚

‚

The first region (Figure 5a), a parcel of intact forest on the border of the states of Zacatecas and
Aguascalientes, is representative of poor model performance with low potential for spurious
change assessments and was selected for its low CI (0.170) and low R2 (0.290). In the land cover
assessments for both the 2001 and 2007, 100% of the pixels were classified as forest, and the
average EVI value over all years was 0.283. Examination of the EVI profiles for 2001 and 2007
indicated similar pattern variability in values over each year.
The second region (Figure 5b), a stretch of irrigated agriculture in the state of Nayarit, at the
very end of the Rio Grande de Santiago delta, is representative of poor model performance with
high potential for spurious change assessments and was selected for its high CI (0.932) and low
R2 (0.336). In both 2001 and 2007 land cover assessments, 100% of the pixels were classified as
agriculture, with an average EVI value of 0.315 for this period. Examination of the EVI profiles for
2001 and 2007 indicated increased vegetative vigor during the growing seasons of July, August,
and September in 2007 compared to 2001.
The third region (Figure 5c), a stretch of seasonal agriculture in the northern portion of the state
of Michoacán near the border of Guanajuato, is representative of high model performance with
low potential for spurious change assessments and was selected for its low CI (0.331) and high
R2 (0.726). In the 2001 land cover map, 49% of the pixels were classified as cropland, and 51%
as mosaic cover. In the 2007 assessment, mosaic cover increased to 72% of the classified pixels,
cropland decreased to 25%, and forest and shrub covered 1% and 2% of the pixels in the parcel,
respectively. The average EVI value for this sample was 0.291 for this period. EVI profiles for 2001
and 2007 demonstrated a very similar pattern for both years.
The fourth region (Figure 5d) is a heterogeneous plot of shrub, forest, and mosaic cover in the
northern region of the state of Jalisco, in an area along a small river basin, where logging is
prevalent. This region demonstrates high potential for spurious change assessments and had
high model performance, as indicated by its high CI (0.955) and high R2 (0.740). In the 2001
land cover map, the parcel was comprised of mosaic (57% of classified pixels), forest (27%),

Remote Sens. 2016, 8, 478

12 of 18

shrub (11%), crop (3%), and water (2%) covers. In 2007, the land cover types included shrub
(58% of classified pixels), mosaic (34%) and forest (8%). The average EVI value from 2001 to 2007
was 0.304. Examination of the average monthly EVI values revealed that values for 2007 were
lower in July than in 2001, but higher in August and September.
4. Discussion
4.1. Variability of Dependent and Independent across Land Cover Types, 2001–2007
Across the Lerma-Chapala-Santiago watershed, the average patterns of EVI, temperature, and
precipitation exhibit noticeable variability through time. Temperature was the most constant variable,
ranging less than 3 degrees C among all years, with the lowest average temperatures (28.5 degrees C)
in 2004 and the highest (31.4 degrees C) in 2001 and 2005. Average EVI increased steadily from
0.255 in 2001 to 0.284 in 2004 before dipping to 0.260 in 2005 and then rising again to 0.276 in 2007.
The precipitation record demonstrates the greatest average differences among years, ranging from a
high in 2004 of 2.58 mm/day to decrease of 70% the following year (0.77 mm/day). Through these
average records, the basic relationship among the independent variables is clear: high precipitation
and low temperatures are associated with high average EVI values.
The relationship among EVI, temperature, and precipitation across each persistent land cover
type related to the typical composition and configuration of the thematic class. Forest, at an average
elevation of 2360.2 m, had the lowest average temperature (25.5 degrees C) and highest average EVI.
Average EVI decreased with the stature of vegetation among natural land covers, based upon the
categorical definitions, and the low average precipitation in the Grass and Sparse classes may indicate
that water stress contributed to the composition of these land covers.
Changing land covers exhibited slightly lower average EVI than persistent land covers, with an
average EVI value of 0.260 for any region that demonstrated a categorical transition and 0.281 for any
persistent land cover. While the high EVI values among persistent land covers may be ascribable to the
intact forest cover and consistent agricultural cultivation from 2001 to 2007, the changing land covers
may also demonstrate high EVI because most transitions in this landscape involved an agricultural
land cover class [1].
The relationship between temperature and EVI is not unidirectional [32], as land cover has a
demonstrable negative effect upon surface temperature [67], but increased temperature may indicate
stresses limiting vegetation [68]. Dense green vegetation, which is characteristic of high average EVI
values, has a lower surface temperature than unvegetated surfaces [32,69]. However, the instantaneous
temperature of the surface of the Earth is a function of the reflected and absorbed solar radiation
and may be influenced by many factors besides land cover type, such as soil moisture, geology, and
topography [63,67]. The positive relationship between temperature and vegetative vigor demonstrated
by others in high latitude and high elevation regions is directly attributable to the limited solar energy
upon vegetation growing patterns [6,11]. In semi-arid regions, such as the Lerma-Chapala-Santiago
watershed, however, the negative relationship between vegetation indices and temperature likely
indicates that water, rather than solar energy, is the greatest limitation on vegetation growth and that
higher temperatures exacerbate dry conditions [7,68].
4.2. The Modeled Relationship of Temperature, Precipitation, and Elevation on EVI
The modeled relationship of EVI as a function of temperature, precipitation and elevation for
each individual year demonstrated higher coefficients of determination than the seven-year series.
Annual models with the highest R2 values were in 2002, 2003, 2004, and 2006 all occurred during
years recognized for El Niño activity [70], typically associated with lower temperature and higher
precipitation. In years recognized for La Niña, which generally brings lower precipitation [70],
72.5% (in 2001) and 68.1% (in 2007) of the variability in EVI values was attributable to climate variables.
The year in which the model performance was lowest, in 2005, with 52.5% of the variability in EVI

Remote Sens. 2016, 8, 478

13 of 18

explained, occurred during a year in which average precipitation indicated drought conditions across
Mexico and in many other regions of the world [39,71].
The average model performance by land cover revealed patterns relating to the dynamics of
vegetation inherent to each cover type in this region. For example, while persistent forest cover had the
highest average EVI during the seven-year period, the model performance was lowest among natural
land covers. Because evergreen pine forests dominate the forest types of this region, the vegetative vigor
measured by a vegetation index demonstrates less interannual variability than broadleaf vegetation.
Natural cover types with the highest R2 between EVI and the independent variables included shrub
(0.605), grass (0.558) and sparse (0.567), indicating that these covers demonstrate a direct and immediate
relationship to changes in precipitation and temperature.
The landscape dynamics of change and persistence of land cover affected the average EVI values.
Regions experiencing a change of categorical land cover between 2001 and 2007 had a closer modeled
response of EVI to independent variables (64.5%) than the average of all persistent land covers, with
only 56.9% explained. Similarly, the mosaic class, defined by the heterogeneous composition of natural
and agricultural land uses showed the highest modeled relationship among all land covers. The
mosaic class comprised more than half of the observed land cover transitions during this time period,
and such transitions often replaced natural land covers with the precipitation-dependent agricultural
activity [27].
4.3. Land Cover Confusion and Natural Variability in Selected Sample Regions
The four sample regions identified through comparison of the Confusion Index and model
performance reveal indicators that can help users of change assessments identify and evaluate the
potential for real and spurious change on the landscape. Previous research identified the prevalence of
change in this landscape using the same data and noted cases where classification uncertainty and
potential error in one or both maps may lead to spurious change assessments [1].
In the forested sample (Figure 5a), a low average CI demonstrates that the signature of this
parcel of land was consistent with the calibration data with which it was classified in the land cover
assessments. As the persistent forest class overall demonstrated low coefficients of determination in
the model, this is consistent with the identification of persistence in this region. These two indices
together indicate that there is little likelihood of change in this parcel.
In the intensive agricultural example (Figure 5b), model performance was relatively low and
confusion was high. This might ordinarily appear to be an indicator of spurious persistence; however
the high CI demonstrates that, in both 2001 and 2007 land cover assessments, this parcel was not very
typical of the range of variability of its thematic category. The legend scheme used for this analysis did
not discriminate between intensive and seasonal agricultural patterns. Because this region is irrigated,
the regression model performance did not demonstrate the close relationship to precipitation and
temperature seen across the rest of the study area. Hence, though this parcel may not be typical of all
pixels in the thematic class, natural variability did not interfere with the change assessment.
Conversely, the change assessment of the seasonal agricultural plot (Figure 5c), with both a low
CI and a high R2 , indicates that natural variability may have a role in the identification of change in
this region. Because seasonal agricultural patterns intermittently leave parcels fallow, a change in the
condition may be perceived as a change in category. The sub-pixel endmembers that comprise the
mosaic class may share similar appearance to both crop or other natural covers, and this class has
been recognized as being very difficult to map accurately [72]. Close examination of the CI sample of
Figure 5c reveals highest confusion in the individual pixels where change was indicated. Though the
land cover assessment indicated change in this region with a low CI value, the high R2 indicates that
the pattern of variability closely follows climatic variables, indicating that the detection of change may
be due to a change in condition rather than a change in category.
Finally, the heterogeneous parcel of mosaic, shrub, and forest cover (Figure 5d), with a very high
CI and high R2 indicates that the while change apparent in the comparison of the 2001 and 2007 land

Remote Sens. 2016, 8, 478

14 of 18

cover maps may not be accurate, something is happening in the region that deviates from normal
patterns of activity. Because both human and natural land covers can show a close relationship to
climate variables, the high R2 alone does not indicate that this region is or is not changing. However,
the high CI of the change assessment, coupled with the increase in EVI in 2007 over 2001 indicates that
the pattern of activity in this region is not typical of other examples of similar classes or other regions
where the model performs well. Indeed this region has undergone tremendous human modification,
including legal and illegal logging as is prevalent across the region [73]. Agriculture is present in
the region, but at an elevation of 917 m over rugged terrain, it is highly dependent on precipitation.
In comparison to the forested example (Figure 5a) which had a low model R2 , the high CI coupled with
high model performance indicates that forests in this parcel maybe be experiencing discrete change
that could be hidden by the natural variability.
5. Conclusions
The relationship of an index of vegetative vigor, such as EVI, to climate variables can reveal
important patterns of variability that can augment land change analyses. The overall performance of
the multiple regression model revealed that an average of 59.1% of the variability in EVI was explained
by concurrent variability in temperature and precipitation, with values as high as 84% explained
variance in some regions. Areas in transition represented a closer link between EVI and the climate
variables than persistent land cover classes, indicating that high variability in vegetation may be
closely linked to land cover change. Further, areas modified by human action for agricultural patterns
(including both seasonal agriculture and the mosaic class) generally demonstrate a closer relationship
between EVI and climate variables than non-agricultural areas. As exemplified by the selected sample
parcels, ultimately, the range of natural variability can both enhance and diminish the perception of
land cover change across a heterogeneous and dynamic landscape, even among confident assessments
of land cover change. Hence, it is essential to consider the vegetative patterns of both natural and
anthropogenic classes in the context of land change assessments in dynamic regions to avoid both
errors of omission and commission.
Acknowledgments: This work was conducted with the support from a NASA Earth System Science Fellowship
(#NNX06AH73H), the Geller Foundation of the Marsh Institute, and the Pruser/Holzhauer Graduate Research
Fund. MODIS and TRMM data were provided by the NASA Earth Observing System Clearinghouse at
http://wist.echo.nasa.gov. SRTM data were acquired from the Consultative Group on International Agricultural
Research Consortium of Spatial Information (CGIAR-CSI) at http://srtm.csi.cgiar.org.
Author Contributions: Zachary Christman prepared data, conducted analysis, interpreted results, and composed
the manuscript and figures, with the advice and guidance of John Rogan, B. L. Turner, II, and J. Ronald Eastman.
Conflicts of Interest: The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:
EVI
MODIS
TRMM
LST
SRTM
AVHRR
NDVI
IGBP-DIS
DEM

Enhanced Vegetation Index
Moderate Resolution Imaging Spectroradiometer
Tropical Rainfall Measuring Mission
Land Surface Temperature
Shuttle Radar Topography Mission
Advanced Very High Resolution Radiometer
Normalized Difference Vegetation Index
International Geosphere-Biosphere Program Data and Information System
Digital Elevation Model

References
1.

Christman, Z.; Rogan, J.; Eastman, J.R.; Turner, B.L. Quantifying uncertainty and confusion in land change
analyses: A case study from central Mexico using MODIS data. GISci. Remote Sens. 2015, 52, 543–570.
[CrossRef]

Remote Sens. 2016, 8, 478

2.
3.
4.
5.
6.
7.
8.

9.
10.

11.

12.
13.
14.
15.
16.
17.
18.
19.
20.

21.
22.
23.

24.

15 of 18

Christman, Z.J. Land Change in Central Mexico: Landscape Heterogeneity, Natural Variability, and Classification
Uncertainty; Clark University: Worcester, MA, USA, 2010.
Escamilla, M.; Kurtycz, A. Social participation in the Lerma-Santiago basin: Water and social life project.
Int. J. Water Resour. Dev. 1995, 11, 457–466. [CrossRef]
Cotler, H. La cuenca Lerma-Chapala: Algunas ideas para un antiguo problema. Gac. Ecol. 2004, 71, 5–10.
Cihlar, J.; St.-Laurent, L.; Dyer, J.A. Relation between the normalized difference vegetation index and
ecological variables. Remote Sens. Environ. 1991, 35, 279–298. [CrossRef]
Schultz, P.A.; Halpert, M.S. Global analysis of the relationships among a vegetation index, precipitation and
land-surface temperature. Int. J. Remote Sens. 1995, 16, 2755–2777. [CrossRef]
Braswell, B.H.; Schimel, D.S.; Linder, E.; Moore, B., III. The Response of Global Terrestrial Ecosystems to
Interannual Temperature Variability. Science 1997, 278, 870–873. [CrossRef]
Rian, S.; Xue, Y.; MacDonald, G.; Touré, M.; Yu, Y.; De Sales, F.; Levine, P.; Doumbia, S.; Taylor, C. Analysis
of climate and vegetation characteristics along the savanna-desert ecotone in mali using modis data.
GISci. Remote Sens. 2009, 46, 424–450. [CrossRef]
Davenport, M.L.; Nicholson, S.E. On the relation between rainfall and the Normalized Difference Vegetation
Index for diverse vegetation types in East Africa. Int. J. Remote Sens. 1993, 14, 2369–2389. [CrossRef]
Gomez-Mendoza, L.; Galicia, L.; Cuevas-Fernandez, M.L.; Magana, V.; Gomez, G.; Palacio-Prieto, J.L.
Assessing onset and length of greening period in six vegetation types in Oaxaca, Mexico, using
NDVI-precipitation relationships. Int. J. Biometeorol. 2008, 52, 511–520. [CrossRef] [PubMed]
Guo, W.Q.; Yang, T.B.; Dai, J.G.; Shi, L.; Lu, Z.Y. Vegetation cover changes and their relationship to climate
variation in the source region of the Yellow River, China, 1990–2000. Int. J. Remote Sens. 2008, 29, 2085–2103.
[CrossRef]
Bradley, B.A.; Mustard, J. Identifying land cover variability distinct from land cover changes: Cheatgrass in
the Great Basin. Remote Sens. Environ. 2005, 94, 204–213. [CrossRef]
Jordan, C.C.F. Derivation of Leaf-Area Index from Quality of Light on the Forest Floor. Ecology 1969, 50,
663–666. [CrossRef]
Sellers, P.J. Relations between canopy reflectance, photosynthesis and transpiration: Links between optics,
biophysics and canopy architecture. Adv. Space Res. 1987, 7, 27–44. [CrossRef]
Tucker, C.J. Red and photographic infrared linear combinations for monitoring vegetation. Remote
Sens. Environ. 1979, 8, 127–150. [CrossRef]
Justice, C.O.; Townshend, J.R.G.; Holben, B.N.; Tucker, C.J. Analysis of the phenology of global vegetation
using meteorological satellite data. Int. J. Remote Sens. 1985, 6, 1271–1318. [CrossRef]
De Beurs, K.M.; Henebry, G.M. Land surface phenology and temperature variation in the International
Geosphere-Biosphere Program high-latitude transects. Glob. Chang. Biol. 2005, 11, 779–790. [CrossRef]
Zhang, X.; Friedl, M.; Schaaf, C.; Strahler, A. Monitoring vegetation phenology using MODIS. Remote
Sens. Environ. 2003, 84, 471–475. [CrossRef]
White, M.A.; Nemani, R.R. Real-Time monitoring and short-term forecasting of land surface phenology.
Remote Sens. Environ. 2006, 104, 43–49. [CrossRef]
Verbesselt, J.; Hyndman, R.; Zeileis, A.; Culvenor, D. Phenological change detection while accounting
for abrupt and gradual trends in satellite image time series. Remote Sens. Environ. 2010, 114, 2970–2980.
[CrossRef]
Townshend, J.R.G.; Justice, C.; Li, W.; Gurney, C.; McManus, J. Global land cover classification by remote
sensing: Present capabilities and future possibilities. Remote Sens. Environ. 1991, 35, 243–255. [CrossRef]
Huete, A.R.; Liu, H.Q.; Batchily, K.; van Leeuwen, W. A comparison of vegetation indices over a global set of
TM images for EOS-MODIS. Remote Sens. Environ. 1997, 59, 440–451. [CrossRef]
Huete, A.R.; Didan, K.; Miura, T.; Rodriguez, E.P.; Gao, X.; Ferreira, L.G. Overview of the radiometric
and biophysical performance of the MODIS vegetation indices. Remote Sens. Environ. 2002, 83, 195–213.
[CrossRef]
Glenn, E.P.P.; Huete, A.R.; Nagler, P.L.L.; Nelson, S.G.G. Relationship between remotely-sensed vegetation
indices, canopy attributes and plant physiological processes: What vegetation indices can and cannot tell us
about the landscape. Sensors 2008, 8, 2136–2160. [CrossRef]

Remote Sens. 2016, 8, 478

25.

26.

27.
28.
29.
30.

31.
32.
33.

34.
35.

36.
37.
38.

39.

40.
41.
42.
43.

44.
45.
46.

16 of 18

Silveira, E.M.D.; de Carvalho, L.M.T.; Acerbi-Junior, F.W.; de Mello, J.M. The assessment of vegetation
seasonal dynamics using multitemporal NDVI and EVI images derived from MODIS. Cerne 2008, 14,
177–184.
Fensholt, R.; Sandholt, I.; Stisen, S. Evaluating MODIS, MERIS, and VEGETATION - Vegetation indices
using in situ measurements in a semiarid environment. IEEE Trans. Geosci. Remote Sens. 2006, 44, 1774–1786.
[CrossRef]
Liu, X.; Kafatos, M. Land-Cover mixing and spectral vegetation indices. Int. J. Remote Sens. 2005, 26,
3321–3327. [CrossRef]
Du Plessis, W.P. Linear regression relationships between NDVI, vegetation and rainfall in Etosha National
Park, Namibia. J. Arid Environ. 1999, 42, 235–260. [CrossRef]
Gao, Z.-Q.; Dennis, O. The temporal and spatial relationship between NDVI and climatological parameters
in Colorado. J. Geogr. Sci. 2001, 11, 411–419.
Kawabata, A.; Ichii, K.; Yamaguchi, Y. Global monitoring of interannual changes in vegetation activities
using NDVI and its relationships to temperature and precipitation. Int. J. Remote Sens. 2001, 22, 1377–1382.
[CrossRef]
Wang, J.; Price, K.P.; Rich, P.M. Spatial patterns of NDVI in response to precipitation and temperature in the
central Great Plains. Int. J. Remote Sens. 2001, 22, 3827–3844. [CrossRef]
Lambin, E.F.; Ehrlich, D. The surface temperature-vegetation index space for land cover and land-cover
change analysis. Int. J. Remote Sens. 1996, 17, 463–487. [CrossRef]
De Beurs, K.M.; Henebry, G.M.; Debeurs, K. Land surface phenology, climatic variation, and institutional
change: Analyzing agricultural land cover change in Kazakhstan. Remote Sens. Environ. 2004, 89, 497–509.
[CrossRef]
West, R.C.; Augelli, J.P. Physical patterns of Middle America. In Middle America, Its Lands and Peoples; Prentice
Hall: Upper Saddle River, NJ, USA, 1976; pp. 22–57.
De Anda, J.; Quinones-Cisneros, S.E.; French, R.H.; Guzmán, M.; de Anda, J.; Quiñones-Cisneros, S.E.;
Guzman, M. Hydrologic Balance of Lake Chapala (Mexico). J. Am. Water Resour. Assoc. 1998, 34, 1319–1331.
[CrossRef]
Magaña, V.O.; Vazquez, J.L.; Perez, J.L.; Perez, J.B. Impact of El Niño on Precipitation in Mexico. Geofica Int.
2003, 42, 313–330.
Cavazos, T.; Hastenrath, S. Convection and rainfall over Mexico and their modulation by the Southern
Oscillation. Int. J. Climatol. 1990, 10, 377–386. [CrossRef]
MacDonald, G.M.; Stahle, D.W.; Villanueva Diaz, J.; Beer, N.; Busby, S.J.; Cerano-Paredes, J.; Cole, J.E.;
Cook, E.R.; Endfield, G.; Gutierrez-Garcia, G.; et al. Climate warming and 21st-century drought in
southwestern north America. EOS Trans. Am. Geophys. Union 2008, 89, 82. [CrossRef]
Stahle, D.W.D.; Cook, E.R.; Villanueva Diaz, J.; Palacio, G.; Fye, F.K.; Burnette, D.J.; Griffin, R.D.; Seager, R.;
Heim, R.R., Jr. Early 21st-century drought in Mexico. EOS Trans. Am. Geophys. Union 2009, 90, 89–100.
[CrossRef]
Pérez-Arteaga, A.; Gaston, K.J.K.J.; Kershaw, M. Undesignated sites in Mexico qualifying as wetlands of
international importance. Biol. Conserv. 2002, 107, 47–57. [CrossRef]
Lind, O.T.; Davalos-Lind, L.O. Interaction of water quantity with water quality: The Lake Chapala example.
Hydrobiologia 2002, 467, 159–167. [CrossRef]
Liverman, D. Vulnerability and adaptation to drought in Mexico. Nat. Resour. J. 1999, 39, 99–116.
Luers, A.L.; Lobell, D.B.; Sklar, L.S.; Addams, C.L.; Matson, P.A. A method for quantifying vulnerability,
applied to the agricultural system of the Yaqui Valley, Mexico. Glob. Environ. Chang. 2003, 13, 255–267.
[CrossRef]
Tucker, C.M.; Eakin, H.; Castellanos, E.J. Perceptions of risk and adaptation: Coffee producers, market shocks,
and extreme weather in Central America and Mexico. Glob. Environ. Chang. 2010, 20, 23–32. [CrossRef]
Farjon, A. Biodiversity of Pinus (Pinaceae) in Mexico: Speciation and palaeo-endemism. Bot. J. Linn. Soc.
1996, 121, 365–384. [CrossRef]
Styles, B.T. Genus Pinus: A Mexican purview. In Biological Diversity of Mexico: Origins and Distribution;
Ramamoorthy, T.P., Bye, R., Lot, A., Fa, J., Eds.; Oxford University Press: Oxford, UK, 1993; pp. 397–420.

Remote Sens. 2016, 8, 478

47.

48.
49.

50.
51.
52.
53.
54.
55.
56.

57.
58.

59.
60.
61.

62.

63.
64.

65.
66.
67.
68.

69.

17 of 18

Brower, L.P.; Castilleja, G.; Peralta, A.; Lopez-Garcia, J.; Bojorquez-Tapia, L.; Diaz, S.; Melgarejo, D.; Missrie, M.
Quantitative changes in forest quality in a principal overwintering area of the monarch butterfly in Mexico,
1971–1999. Conserv. Biol. 2002, 16, 346–359. [CrossRef]
Lopez, E.; López, E.; Bocco, G.; Mendoza, M.; Duhau, E. Predicting land-cover and land-use change in the
urban fringe: A case in Morelia city, Mexico. Landsc. Urban Plan. 2001, 55, 271–285. [CrossRef]
Alvarez, R.; Bonifaz, R.; Lunetta, R.S.C.; García, C.; Gómez, G.; Castro, R.; Bernal, A.; Cabrera, A.L.
Multitemporal land-cover classification of Mexico using Landsat MSS imagery. Int. J. Remote Sens. 2003, 24,
2501–2514. [CrossRef]
Von Bertrab, E. Guadalajara’s water crisis and the fate of Lake Chapala: A reflection of poor water
management in Mexico. Environ. Urban 2003, 15, 127. [CrossRef]
Levine, G. The Lerma-Chapala river basin: A case study of water transfer in a closed basin. Paddy
Water Environ. 2007, 5, 247–251. [CrossRef]
Klooster, D. Campesinos and Mexican forest policy during the twentieth century. Lat. Am. Res. Rev. 2003, 38,
94–126. [CrossRef]
Avalos, H.C.; Santander, A.G.P.; Rodríguez, C.; Guadarrama, C.E. Determinación de zonas prioritarias para
la eco-rehabilitación de la cuenca Lerma-Chapala. Gac. Ecol. 2004, 71, 79.
Sanderson, S.E. The Transformation of Mexican Agriculture: International Structure and the Politics of Rural Change;
Princeton University: Princeton, NJ, USA, 1986.
Appendini, K.; Liverman, D. Agricultural policy, climate change and food security in Mexico. Food Policy
1994, 19, 149–164. [CrossRef]
Zahniser, S.; Link, J. Effects of North American Free Trade Agreement on Agriculuture and the Rural Economy;
Agriculture and Trade Report No. (WRS-0201). Economic Research Service, United States Department of
Agriculture: Washington, DC, USA, 2002.
Dalton, R. Saving the agave. Nature 2005, 438, 1070–1071. [CrossRef] [PubMed]
Loveland, T.R.; Reed, B.C.; Brown, J.F.; Ohlen, D.O.; Zhu, Z.; Yang, L.; Merchant, J. Development of a global
land cover characteristics database and IGBP DISCover from 1 km AVHRR data. Int. J. Remote Sens. 2000, 21,
1303–1330. [CrossRef]
Friedl, M.; Strahler, A.; Zhang, X.J. The MODIS land cover product: Multi-attribute mapping of global
vegetation and land cover properties from time series MODIS data. Geosci. Remote 2002, 6, 3199–3201.
Olson, J.S. Global Ecosystems Framework: Definitions; Internal Report; USGS EROS Data Center: Sioux Falls,
ND, USA, 1994.
Wolfe, R.E.; Roy, D.P.; Vermote, E. MODIS land data storage, gridding, and compositing methodology: Level
2 grid. (Moderate Resolution Imaging Spectroradiometer)(Special Issue on EOS AM-1 Platform, Instruments,
and Scientific Data). IEEE Trans. Geosci. Remote Sens. 1998, 36, 1324–1338. [CrossRef]
Huffman, G.J.; Adler, R.F.; Bolvin, D.T.; Gu, G.; Nelkin, E.J.; Bowman, K.P.; Hong, Y.; Stocker, E.F.; Wolff, D.B.
The TRMM Multisatellite Precipitation Analysis (TMPA): Quasi-Global, Multiyear, Combined-Sensor
Precipitation Estimates at Fine Scales. J. Hydrometeorol. 2007, 8, 38–55. [CrossRef]
Lookingbill, T.R.; Urban, D.L. Spatial estimation of air temperature differences for landscape-scale studies in
montane environments. Agric. For. Meteorol. 2003, 114, 141–151. [CrossRef]
Huang, S.L.; Rich, P.M.; Crabtree, R.L.; Potter, C.S.; Fu, P.D. Modeling monthly near-surface air temperature
from solar radiation and lapse rate: Application over complex terrain in Yellowstone National Park.
Phys. Geogr. 2008, 29, 158–178. [CrossRef]
Christman, Z.; Rogan, J. Error propagation in raster data integration: Impacts on landscape composition and
configuration. Photogramm. Eng. Remote Sens. 2012, 78, 617–624. [CrossRef]
Ronald Eastman, J.; Sangermano, F.; Ghimire, B.; Zhu, H.; Chen, H.; Neeti, N.; Cai, Y.; Machado, E.; Crema, S.
Seasonal trend analysis of image time series. Int. J. Remote Sens. 2009, 30, 2721–2726. [CrossRef]
Buyantuyev, A.; Wu, J.G. Urban heat islands and landscape heterogeneity: Linking spatiotemporal variations
in surface temperatures to land-cover and socioeconomic patterns. Landsc. Ecol. 2010, 25, 17–33. [CrossRef]
Karnieli, A.; Agam, N.; Pinker, R.T.T.; Anderson, M.; Imhoff, M.L.L.; Gutman, G.G.G.; Panov, N.; Goldberg, A.
Use of NDVI and land surface temperature for drought assessment: Merits and limitations. J. Clim. 2010, 23,
618–633. [CrossRef]
Friedl, M.; Davis, F. Sources of variation in radiometric surface temperature over a tallgrass prairie.
Remote Sens. Environ. 1994, 48, 1–17. [CrossRef]

Remote Sens. 2016, 8, 478

70.

71.
72.
73.

18 of 18

National Weather Service Climate Prediction Center. Cold & Warm Episodes by Season. Available online:
http://www.cpc.ncep.noaa.gov/products/analysis_monitoring/ensostuff/ensoyears.shtml (accessed on
30 March 2016).
Boyd, R.; Ibarrarán, M.E. Extreme climate events and adaptation: An exploratory analysis of drought in
Mexico. Environ. Dev. Econ. 2009, 14, 371–395. [CrossRef]
Loveland, T.R.; Belward, A.S. The IGBP-DIS global 1 km land cover data set, DISCover: First results. Int. J.
Remote Sens. 1997, 18, 3289–3295. [CrossRef]
Honey-Rosés, J. Disentangling the proximate factors of deforestation: The case of the Monarch butterfly
Biosphere Reserve in Mexico. Land Degrad. Dev. 2009, 20, 22–32. [CrossRef]
© 2016 by the authors; licensee MDPI, Basel, Switzerland. This article is an open access
article distributed under the terms and conditions of the Creative Commons Attribution
(CC-BY) license (http://creativecommons.org/licenses/by/4.0/).