Introduction

Ecologists have long debated the relationship between environmental heterogeneity and biodiversity (e.g., MacArthur and Pianka 1966; Wiens 2002; Seiferling et al. 2014). The consensus describes environmental heterogeneity as being an important driver of biodiversity maintenance and ecosystem health. For rapid assessment in biodiversity monitoring, many studies have, therefore, focussed on modelling relationships between landscape heterogeneity and species diversity (e.g., MacArthur and MacArthur 1961; Tews et al. 2004; Tamme et al. 2010; Zhao et al. 2015). To this end, environmental heterogeneity is often considered equivalent to landscape heterogeneity in practice (Tscharntke et al. 2012). However, such practice inevitably reduces the complex realism of environmental heterogeneity into discrete patches that may not physically or functionally occur in reality (Turner 1989; Cushman et al. 2010; Fahrig et al. 2011). Although such a mosaic approach of discretizing environmental/landscape heterogeneity has been successful, especially in urban and agricultural landscapes, it falls short in natural ecosystems where the classification of these patches discounts important within-patch heterogeneity (McGarigal et al. 2009).

Environmental gradients are considered as an alternative to this mosaic approach, one which arguably better reflects the continuous nature of environmental heterogeneity (Doebeli and Dieckmann 2003; Guisan and Thuiller 2005). Remote sensing offers a cost-effective, systematic and repeatable method of mapping and monitoring environmental heterogeneity as a continuous surface (e.g., González-Megías et al. 2011; Hernández-Stefanoni et al. 2012; Duro et al. 2014). The spectral response of satellite imagery is therefore often used to analyse ecosystem patterns and processes (Gould 2000; Wulder et al. 2004). Variations in this spectral response can originate from corresponding variations in the underlying properties of the physical landscape as well as other biological features (Rocchini et al. 2013). Separating out the different drivers of environmental heterogeneity from this spectral response however remains a challenge (Somers et al. 2011; Shi and Wang 2014).

There is to date no definitive method to quantify environmental heterogeneity, as such, a robust environmental heterogeneity–biodiversity relationship remains elusive (Allouche et al. 2012; Redon et al. 2014). Recent studies further suggest that the relationship itself is non-ubiquitous (Bar-Massada and Wood 2014), varying across scale (Oldeland et al. 2010; Stein et al. 2014), level of ecosystem disturbance (Seiferling et al. 2014), species geographic range (Katayama et al. 2014) and available habitat area (Fahrig 2013). We expect this is due to the complexity of environmental heterogeneity and the contingency of identifying key drivers of this heterogeneity using conventional methods (Johnson 2007). Nevertheless, in the face of increasing concerns of global biodiversity loss (MEA 2005; Hooper et al. 2012) how environmental heterogeneity is defined and measured is a key question for today’s conservation agencies. We define environmental heterogeneity here as the variation of landscape form (the physical rendition of composition, structure and function), as represented by Landsat spectral variation. With the term landscape complexity we refer to the interacting processes that are underlying the observed environmental heterogeneity. In this paper we operationalize landscape complexity by measuring what part of the environmental heterogeneity is not straightforwardly explained by spatial and temporal variation in the biophysical context.

While many studies have sought to develop cost-effective, systematic and repeatable methods of mapping and monitoring biodiversity (e.g., Duro et al. 2007; Reyers and McGeoch 2007; Lengyel et al. 2008; Pettorelli et al. 2014), few have explored the spatial and temporal variability of environmental heterogeneity itself. Using traditional global models, inherent spatial structures are often ignored and important information about how relationships between observed heterogeneity and physical landscape properties might change over space discounted (Guo et al. 2008; Matthews and Yang 2012). For instance, in the Kruger National Park (KNP) in South Africa, we would expect to find highly variable relations between environmental heterogeneity and other spatially explicit drivers, such as geology and soils, topography, climate, herbivory and fire. Geographically weighted regression (GWR) is reportedly able to incorporate these local spatial relationships into a traditional regression framework (Brunsdon et al. 2002; Fotheringham et al. 2002). In contrast to global regression techniques, GWR accounts for spatial non-stationarity by allowing relationships to vary over space. In this way, if the elicited response varies geographically it would suggest different processes are interacting within the study area (Matthews and Yang 2012). We therefore anticipate GWRs application in KNP to enable us to visualise the geographical variation of environmental heterogeneity and identify its key drivers across the park (Oliveira et al. 2014).

In this paper, we use local GWR models to map the relationship between Landsat spectral variance and stable physical landscape properties. We explore how this relationship changes over space and time and examine how in some areas heterogeneity patterns may be driven more noticeably by dynamic ecological processes. We used spectral variation as a proxy for environmental heterogeneity which depicts the variability of a spectral response across different wavelengths or bands of a Landsat satellite image (Short 2005). For stable physical landscape properties we used landscape features that do not change over ~50 years, such as elevation and geology. Based on our findings, we identify the proportion of spectral variability in the landscape, as seen from the multispectral Landsat satellite, explained by stable landscape properties. Thereafter, we examine the sensitivity of this relationship to changes in season and rainfall and explore how mapped model outcomes change as a result. We propose that, the proportion of model disagreement (1−R2), represents the level of landscape complexity, i.e. the influence of dynamic landscape processes and stochastic disturbance events, not fully captured by KNPs underlying physical landscape template. We investigate how landscape complexity has changed spatially in KNP over time and highlight areas where the change in this complexity has been consistent, potentially signalling rapid changes in underlying dynamic ecological processes that could strongly affect biodiversity. To this end, we test the degree to which landscape complexity can explain local plant species richness patterns and provide insight into its application for protected area managers.

Methods

Study area

KNP has considerable biophysical diversity and a long conservation history (du Toit et al. 2003). It is one of the largest protected areas (PAs) in the world (~2 m ha), situated in the north-eastern corner of South Africa (Fig. 1). The park is dominated by gently undulating topography (150–840 m a.s.l) underlying granite gneiss, schists, amphibolites, basalt and gabbros (Schutte 1986). Mountainous areas occur in the east, along the border of Mozambique (Lebombo Mountains), in the south-west (Malelane Mountains) and in the north-west (Soutspansberg Mountains) (Schutte 1986). A basalt-granite east–west division is clearly visible with the more fertile basalts in the east and less fertile granites in the west (Munyati and Ratshibvumo 2010). Climate is a major ecosystem driver (Pickett et al. 2003; Venter et al. 2008) with rainfall occurring in decadal wet and dry cycles across KNP. The long-term annual mean varies from 350 mm in winter to 950 mm in summer with a slight dry–wet rainfall gradient occurring from north to south-west (Gertenbach 1980). Average temperatures range from 26.4 °C in summer (December–March) to 17.8 °C in winter (June–August) (Zambatis 2006). KNP falls within South Africa’s dominant savanna biome (Low and Rebelo 1996), an inherently heterogeneous ecosystem driven by complex spatial interactions between rainfall, soil, vegetation patterns and dynamic processes such as herbivory, fire and floods (Groen 2007).

Fig. 1
figure 1

Kruger National Park, situated in the north-eastern corner of South Africa between latitudes 22°19′40″S–25°31′44″S and longitudes 30°53′18″E–32°01′59″ within the country’s dominant Savanna Biome, overlaying a gently undulating topography

Full size image

Data analyses

All analyses were carried out in R version 3.0.2 (R Core Team 2013), RStudio version 0.98.978 (RStudio 2013) and GRASS GIS version 7.1.svn (GRASS 2014) in a step-wise manner: (1) Landsat spectral variation, (2) physical landscape variation, (3) GWR models and our interpretation of model fit in terms of landscape complexity, (4) Landscape complexity’s relationship to plant species richness.

Landsat spectral variation

Landsat imagery were available for the Skukuza region (path 168—row 077, WRS2) from different sensors (MSS, TM, ETM+, OLI) since 1972. Six images, representing late season winter conditions (July or August months) and summer conditions (March or April months), were selected for years signifying long-term mean, below and above average rainfall periods. Representative years were selected using a three-year rolling mean of daily rainfall records from the Skukuza weather station and associated availability of cloud-free images (Fig. S1). Final image dates represent winter and summer ‘windows’ into low (1991-07-30; 1993-04-14), average (1984-08-27; 1987-03-13) and high (1998-08-18; 2000-04-09) rainfall conditions. Inherent sources of error were dealt with as follows: digital numbers were converted into surface reflectance units using the US Geological Survey’s (USGS) on demand interface for the Earth Resources Observation and Science’s (EROS) Centre Science Processing Architecture (ESPA Ordering Interface 2013); each band (excluding band 6) was geometrically and radiometrically calibrated to the standard terrain correction (1T) level (Irish 2000) with a UTM WSG84 36S projection using GRASS (2014).

Before calculating spectral variation, a correlogram (Wright 2015) and local Moran’s I measure of spatial autocorrelation (Hijmans 2015) were calculated for bands 1–5 and 7, revealing a non-stationary covariance structure, typical of remote sensing data (Wulder and Boots 2000; Propastin 2009). That is, bands were found to be significantly collinear (Fig. S2) and spatially autocorrelated (Table S1). We removed the first source of error (inter-band collinearity) by transforming individual bands into principle components (PC) using the i.pca function in GRASS (2014) interfaced through R (R Core Team 2013). Spatial autocorrelation (i.e. intra-band collinearity) was addressed through the use of GWR (discussed later).

Once individual bands were transformed into PCs, the resulting eigenvalues (or loadings; summarised in Table S2) explained the proportion of variance accounted for by each PC across the different years. For example, a high PC1 loading would suggest a large percentage of the variation in the landscape can be measured using only the first principle axis (Ringnér 2008). Conversely, a low PC1 loading would suggest one axis rotation is not enough to account for all the variability in the landscape and therefore the structure of the data, and in our case the landscape, would be more complex. Exploratory results indicate the proportion of variance accounted for by PC1, for example, is generally higher in winter and lower rainfall periods compared to summer and higher rainfall periods (Fig. S3). This suggests season and rainfall are potentially important drivers of landscape complexity. To better understand environmental heterogeneity (Rocchini and Neteler 2012), we further calculated the textural variance, entropy and uniformity (or angular second moment (ASM)) for each PC (1–6) using r.texture within a three by three pixel moving-window neighbourhood, as well as Shannon’s, Simpson’s and Rényi’s Entropy diversity indices and Pielou’s Evenness index using r.diversity [See Rocchini et al. (2013) and the GRASS (2014) reference manual for details about index formulas].

Physical landscape variation

Environmental heterogeneity depicted by spectral variation is then regressed by the variability of stable physical landscape properties. These are underlying properties of the landscape template which do not change over ~50 years, namely elevation, slope, aspect, flow direction, watershed area, potential surface wetness index and soil form, depth and clay content (Fig. 2).

Fig. 2
figure 2

Stable physical landscape elements that do not change over a 50 year period, forming Kruger National Park’s physical landscape template: a elevation, b slope, c aspect, d watershed area, e potential surface wetness index, f soil form, g soil depth and h soil clay content. Profile graphics in the margins illustrate mean latitudinal and longitudinal values (Perpiñán and Hijmans 2014)

Full size image

KNPs slope and aspect were calculated from a 5 m digital elevation model (DEM) (Van Niekerk 2012) using r.slope.aspect in GRASS (2014); flow direction, watershed area (sink) and a surface wetness index (or topographic convergence index (TCI)) using r.terraflow (GRASS 2014). Soil form, depth and clay content were extracted from the Mpumalanga Province Natural Resources dataset (Wessels et al. 2001). We selected an uncorrelated subset of explanatory variables using the variance inflation factor (VIF), which excludes highly correlated variables through a stepwise procedure (Naimi 2015). Flow direction, which was negatively correlated with aspect (−0.63) and soil clay content, which was positively correlated with soil form (0.67), both had higher VIF values and were therefore removed along with watershed area (VIF = 3.1). After removing these variables, final VIF scores were satisfactorily all below 1.5 (Fig. S4). We continued with our analysis using elevation, aspect, slope, TCI, soil form and soil depth as our explanatory variables. As with spectral variation, we express their variability in the landscape in terms of both textural features measured as variance, entropy and uniformity (r.texture) as well as the same diversity indices (r.diversity) of properties (Rocchini et al. 2013), within a three by three moving window area (GRASS 2014).

Geographically weighted regression (GWR) and landscape complexity

The relationships between the resulting measures of variance for spectral and physical landscape properties were estimated using GWR (Gollini et al. 2015) for different seasons and rainfall conditions. We included season and rainfall because they could potentially affect vegetation structure, the intensity of disturbance (e.g., fires) and the distribution of large fauna (Chirima et al. 2012; Smith et al. 2013). Using the R package GWmodel, we identified an optimal bandwidth for each model based on the Akaike Information Criterion (AIC) with an adaptive bisquare bandwidth setting (Gollini et al. 2015) (Table S3). We compared model fit of the different variance measures using AIC and selected the ‘best’ measure to explore relationships further. Thereafter, we examined how GWR model results vary across winter and summer months of representative low (1991–1993), average (1984–1987) and high (1998–2000) rainfall periods using a multiple comparison test after Kruskal–Wallis (Giraudoux 2015) as well as an analysis of variance model (ANOVA). We then calculated contrasts for factor interactions to explore how seasonal contrasts of GWR coefficients differ between rainfall groups (de Rosario-Martinez 2015). Resulting local adjusted coefficients of determination (R2) were mapped to highlight the spatial variability of model performance against season and rainfall. Spatial non-stationarity was tested using Leung’s F3 statistic (Leung et al. 2000). Regressions were run on a sample (n = 2586) of the original raster data.

We qualified landscape complexity as the level of model disagreement (1−R2), representing the relationship of spectral response to dynamic processes not fully captured by the underlying stable landscape template alone. We examined how landscape complexity changed over time by detecting the trend in 1−R2 values from 1984 until 2000 across the surface of KNP. We expect a change in landscape complexity to be indicative of changes in the driving ecological factors behind environmental heterogeneity. Results are summarised as surface trend maps indicating areas in the KNP where the degree of change in landscape complexity fluctuates, and thus possibly biodiversity, with changing seasonal and rainfall conditions.

Landscape complexity and plant species richness

The proportion of total variation in spectral response explained by physical landscape properties is captured by the R2 from GWR and is a measure of model agreement. The remaining, unexplained proportion (1−R2) therefore represents spectral variation that cannot be explained by physical landscape properties alone. There are many dynamic landscape properties that could help explain this remaining variation in the landscape, for example fire, vegetation dynamics, herbivore distribution and a human footprint. However, detailed records of these properties are rarely available. As an alternative, we interpret 1−R2 as a measure of landscape complexity, distinguishing the level of influence of dynamic landscape processes and stochastic disturbance events, from the underlying physical landscape template.

We tested this theory by examining the degree to which landscape complexity explained local patterns of plant species richness. Woody plant species data were obtained from the historical surveys of Venter (1990), recently described by Kiker et al. (2014). These data contain detailed surveys of woody vegetation cover and composition subset to our study area (n = 692 sites, totalling 115 species). A species accumulation curve (SAC) was computed (Oksanen et al. 2015) using the random method to find mean SAC and the number of species for all sample sites in our study area. Relationships between resulting species richness per site and landscape complexity were assessed using, again, GWR. We summarised GWR results and examined how parameter estimates vary with season (winter and summer) and rainfall (low, average and high conditions).

Results

Spectral and physical landscape variation

Box-and-whisker diagrams illustrate the shape of variation in spectral response of Landsat PCs across seasons and a rainfall gradient (Fig. 3). In general, dispersion of PC values tends to increase with increasing rainfall in winter months but decreases as rainfall increases in summer months (Fig. 3a). Textural measures of randomness (entropy, Fig. 3c) and its converse, uniformity (Fig. 3d), showed similar seasonal patterns, i.e. winter entropy increased while summer entropy decreased and winter uniformity decreased while summer uniformity increased as rainfall increased. Diversity clearly increased as rainfall increased across both winter and summer months (Kruskal–Wallis χ2 = 13792.39, df = 5, p-value < 0.0001; Fig. 3e–g). Variability of physical landscape properties (elevation, slope, aspect, flow direction, watershed area, potential surface wetness index and soil form, depth and clay content) are unchanged by year or season.

Fig. 3
figure 3

Boxplots assessing the location, dispersion, and symmetry or skewness of spectral variation, as measured by different indices, across different seasons and rainfall conditions: a Raw principle components (PC) of Landsat bands 1–5 and 7 (Kruskal–Wallis χ2 (K–W χ2) = 30081.27, df = 5, p < 0.0001); b PC textural variance (K–W χ2 = 11420.06, df = 5, p < 0.0001; c PC textural entropy (K–W χ2 = 12332.45, df = 5, p < 0.0001; d PC textural uniformity (K–W χ2 = 12189.02, df = 5, p < 0.0001; e PC Shannon’s diversity (K–W χ2 = 10276.83, df = 5, p < 0.0001); f PC Simpson’s diversity (K–W χ2 = 13792.39, df = 5, p < 0.0001); g PC Rényi’s diversity (K–W χ2 = 11715.88, df = 5, p < 0.0001; h PC Pielou’s evenness (K–W χ2 = 1637.187, df = 5, p < 0.0001). Brackets indicate differences which are not significant according to the Kruskal–Wallis rank sum test (Giraudoux 2015). All other differences are significant

Full size image

Geographically weighted regression (GWR) and landscape complexity

Models with raw PC values representing spectral variation (response variables) and raw physical landscape properties (explanatory variables) were consistently better able to balance model fit and complexity than other indices, as indicated by the notably lower AIC scores (Fig. S5). Therefore, we only examined the local relationships between the linear combination of spectral PC values and uncorrelated stable physical landscape properties.

GWR results show that relationships between spectral variation and KNPs physical landscape template changed with season and rainfall and were spatially diverse (Table S3). Leung et al.’s (2000) F3 test for spatial non-stationarity shows elevation, surface wetness and soil form estimates vary significantly over the region for all years (Table 1). Whereas aspect, slope and soil depth appear constant in some years but vary significantly in others (Table 1). The range of R2 values was wide (0.1–0.7), varying considerably across the landscape and over time (Fig. 4). The proportion of spectral variance captured by physical landscape properties, as described by R2, also varied within and between years (Fig. 4). A multiple comparison test after Kruskal–Wallis (Giraudoux 2015) indicated season and rainfall class both had a significant effect on R2 values (Kruskal–Wallis χ2 = 11951.46, df = 5, p < 0.0001).

Table 1 Significance of non-stationarity in the physical landscape variable’s coefficient estimates the GWRs after Leung et al. (2000) (see Table S3 in supplementary material for full results)
Full size table
Fig. 4
figure 4

Geographically weighted regression (GWR) R2 values for winter and summer months of representative low (1991–1993), average (1984–1987) and high (1998–2000) rainfall periods. Maps show the spatial heterogeneity in the proportion of spectral variance accounted for by physical landscape properties with season and rainfall (Kruskal–Wallis χ2 = 11951.46, df = 5, p < 2.2e-16). ANOVA results show model fit (R2) generally increased from low to high rainfall (1–3; β = 0.06, t(36,420) = 34.184, p < 0.001) and from winter to summer (a to b; β = 0.05, t(36,420) = 28.573, p < 0.001). R2 values were significantly lower in higher rainfall summer months (b3) compared to lower rainfall winter months (a1) (β = −0.17, t(36,420) = −73.479, p < 0.001). A contrast interaction test (de Rosario-Martinez 2015) showed R2 seasonal contrasts differed significantly between rainfall groups: i.e. summer R2 low (b1) to high rainfall (b3) contrasts were 0.17 less than those in winter months (a1 and a3) (β = −0.174373, df = 1, SS = 46.149, F = 5399.2, p < 0.0001). Similarly low rainfall R2 winter (a1) to summer (b1) contrasts were 0.26 less than those for high rainfall periods (a3 and b3) (β = −0.259394, df = 2, SS = 106.16, F = 6210.2, p < 0.0001)

Full size image

On the surface, GWR results show model fit (R2) generally increased from low to high rainfall (β = 0.06, t(36,420) = 34.184, p < 0.001) and from winter to summer (β = 0.05, t(36,420) = 28.573, p < 0.001). However, when adding an interaction effect between season and rainfall, this result was reversed for summer months. That is, R2 values were significantly lower in higher rainfall summer months compared to lower rainfall winter months (β = −0.17, t(36,420) = −73.479, p < 0.001). A contrast interaction test (de Rosario-Martinez 2015) confirmed R2 seasonal contrasts differed significantly between rainfall groups: i.e. summer R2 low to high rainfall contrasts were 0.17 less than those in winter months (β = −0.174373, df = 1, SS = 46.149, F = 5399.2, p < 0.0001). Similarly low rainfall R2 winter to summer contrasts were 0.26 less than those for high rainfall periods (β = −0.259394, df = 2, SS = 106.16, F = 6210.2, p < 0.0001).

The spatial trend surface of landscape complexity from 1984 to 2000 (Fig. 5), revealed areas (in green) where landscape complexity (1−R2) has increased, areas (in red) where it has declined or areas (in yellow) where it has remained relatively unchanged from 1984 to 2000.

Fig. 5
figure 5

Total accumulated difference of landscape complexity (GWR 1−R2) represented as surface trends across low (1991–1993), average (1984–1987) and high (1998–2000) rainfall periods. Shades of green indicate areas that increased in landscape complexity, shades of red decreased in landscape complexity and shades of yellow remain unchanged from 1984 to 2000

Full size image

Landscape complexity and plant species richness

GWR R2 results mapped over the spatial extent of our study area (Fig. 4) illustrate the degree to which model agreement differed spatially across winter and summer months of representative low (1991–1993), average (1984–1987) and high (1998–2000) rainfall periods. We interpret its inverse, 1−R2 (model disagreement), as the level of complexity in the landscape. GWR results show a significant proportion of the variance in plant species richness can be explained by our measure of landscape complexity (R2 values ranged from 0.70 to 0.78) (Table 2). These results showed significant improvement over GWRs of raw physical landscape properties and raw surface reflectance PC values, which only accounted for 62 and 57 % of the variance in plant species richness respectively (Table S4). There was a significant positive correlation between plant species richness and landscape complexity in the years closest to sample collection dates ~1989 (1987: R2 = 0.74, median β = 0.35, F(1691) = 21.36, p < 0.0001; 1993: R2 = 0.78, median β = 6.43, F(1691) = 18.07, p < 0.0001). A map of the residuals of plant species richness for 1993 illustrates the spatial variability of this relationship (Fig. 6). In contrast, a lower but still significant negative correlation between plant species richness and landscape complexity was found in 1991 (R2 = 0.76, median β = −7.46, F(1691) = 4.623, p = 0.0319), which could be explained by the confounding effects of the severe drought KNP experienced in 1991/1992 (Zambatis and Biggs 1995).

Table 2 GWR results of plant species richness modelled as a function of model fit (R2)
Full size table
Fig. 6
figure 6

A map of the residuals (observed—fitted plant species richness) for 1993 illustrating the spatial variability of the relationship between plant species richness with landscape complexity. Results described in Table 2 show an overall significant positive correlation varying spatially between plant species richness and landscape complexity (R2 = 0.78, median β = 2.45, F(1691) = 18.07, p < 0.0001)

Full size image

Discussion

In their meta-analysis, Stein et al. (2014) found environmental heterogeneity to be an important driver of species richness. Remotely sensed spectral heterogeneity is recommended by several authors (Duro et al. 2007; Rocchini et al. 2010; Nagendra et al. 2013; Pettorelli et al. 2014) as a proxy for environmental heterogeneity and the consequent rapid assessment of biodiversity properties. It stands to reason that, a more diverse spectral response will represent a more diverse landscape in that spectral heterogeneity will reflect the associated variation of environmental properties in the landscape. However, we also expect this relationship to be dynamic, changing across different and interactive space-time scales. We demonstrated this using Fotheringham et al.’s (2002) GWR technique with Landsat surface reflectance and stable physical landscape properties. By allowing relationships to vary over space, we were able to account for spatial non-stationarity and visualise the resulting patterns (Brunsdon et al. 1996). Such ability is especially important for ecological studies where the geographic-structure of relationships between predictors and response variables are likely to manifest differently in space. Windle et al. (2010), for example, have shown GWR to be superior to other global methods commonly used in terrestrial ecology, including, generalised linear, additive, and linear mixed models. However, ecologists should also be weary of the limitations of GWR, such as the effects of collinearity highlighted by Finley (2011) and Czarnota et al. (2015), and the difficulties of appropriate bandwidth selection described by Matthews and Yang (2012). Pasher et al. (2013) further describe these in their approach to improve estimates of landscape effects on ecological responses. We dealt with this using principle components analysis, correlograms, VIFs and an adaptive bandwidth selection using AICs.

Results showed that the relationship between spectral heterogeneity and stable physical landscape properties is sensitive to season and rainfall condition. Moreover, we showed that textural measures of entropy increased with rainfall in winter but decreased with rainfall in summer. While, textural measures of uniformity (ASM) also showed an inverse pattern of decreasing uniformity with increasing rainfall in winter and increasing uniformity with increasing rainfall in summer. We suggest these results are representative of both (1) true structural diversity in the landscape and (2) the limitations of remotely sensed Landsat data. In the first instance, we expected structural diversity to increase with rainfall up to a threshold where vegetation cover, for example, would reach an asymptote thereby decreasing structural entropy and increasing structural uniformity, as seen in Fig. 3. However, in the second we also recognise that this outcome may be an effect of what is ‘visible’ to the Landsat’s passive sensor. Under dense and extensive cover conditions, this satellite is less able to detect under-canopy variability in the landscape. Prospective studies may wish to explore the use of active sensors like Lidar in future.

Regionally, spectral diversity increased with increasing rainfall across both winter and summer months. Intuitively these results represent the increase in environmental diversity as water availability becomes less limiting. This is corroborated by our findings that the proportion of satellite surface reflectance variance captured by a single PC axis rotation, for example, was generally higher in winter and during lower rainfall periods as compared to summer and higher rainfall periods. Locally, raw PC values representing spectral variation and raw physical landscape properties were consistently better able to balance model fit and complexity than other textural or diversity measures. This is consistent with the findings of Warren et al. (2014), who found spectral diversity yielded reasonable estimates of plant species richness using a simple Pearson correlation to measure linear relationship strength. Our results add a spatial component which also proved spatial non-stationarity was statistically significant. For ecological studies this highlights the importance of using a geographical approach when analysing environmental data to capture different responses driven by non-stationary processes (Leung et al. 2000; Brunsdon et al. 2002; Matthews and Yang 2012). This was made further evident by the clear differences in the proportion of spectral variation captured by relatively stable physical landscape properties over space and time. Spatially, the spectral variation explained by stable physical properties (R2) varied widely across the landscape. This reiterates that environmental heterogeneity is driven not only by stationary physical landscape properties but also non-stationary or dynamic processes. These dynamic processes are often difficult to isolate but their compound influence on the landscape can be measured as shown by the wide range of R2 values; or its inverse (1−R2) that is landscape complexity.

Over time, we found that increasing summer rainfall reduces the explanatory power of stable physical landscape properties on environmental heterogeneity (as measured by Landsat spectral variation). We postulate that this general reduction in the explanatory power of models fitted to data from summer periods, and periods of high rainfall versus winter and lower rainfall periods, is indicative of dynamic environmental processes not captured by physical landscape properties. These dynamic processes are driven by season and rainfall and include, for example fire, vegetation dynamics, herbivore distribution, and human development. Under higher rainfall conditions vegetation activity, for instance, is increased and herbivore density and distribution patterns will change in response.

We hypothesised that the proportion of spectral variation unexplained by the underlying physical landscape template is representative of the level of complexity in the landscape. Little temporal trend in landscape complexity between seasons (winter-summer) and rainfall (low-average-high) conditions could potentially highlight comparatively stable landscapes (Fig. 5). Changing ecological drivers in the landscape could show a more consistent change in landscape complexity. Geographically these areas tend to coincided with basalt dominated areas in the east (Fig. 5), suggesting KNP basalts are generally becoming more complex than their granitic counterparts in the west. This is consistent with the findings of Colgan et al. (2012) who showed above-ground biomass production on basalts was driven largely by herbivore-fire interactions rather than soil properties. In other words, basaltic areas appear to be driven more by dynamic ecological processes and feedbacks.

We tested this theory against plant species richness data and found a strong, significant relationship between landscape complexity and species richness, with areas presenting negative residuals of species richness potentially associated with higher elevation and granites (Fig. 6). These findings show that indeed processes other than physical landscape properties shape environmental heterogeneity and biodiversity over space and time. In support of this, Proulx et al. (2015) also found a variety of climate-biodiversity relationships. They suggest drivers of biodiversity are built on complex interactions of environmental, within-species and between-species variability (Proulx et al. 2015). We suggest plant species richness in the KNP shows a stronger relationship with landscape complexity than physical landscape properties because of this multiplicity of effects.

However, empirical knowledge of dynamic processes is often not available for protected area managers, and even if accessible, is rarely spatially explicit or temporally continuous. Nevertheless, such knowledge remains central to understanding the functioning of natural systems and their effectual management as protected areas. Our approach provides a starting point by mapping the relative importance of stable physical landscape properties compared to other unknown dynamic processes for environmental heterogeneity. We showed how dynamic processes move across the landscape over time and suggest that biodiversity monitoring programmes be designed to capture this variability. For example, KNP’s annual herbivore counts, done solely in the dry-season, may be missing important changes in distribution patterns driven by seasonal changes in landscape complexity (Martin et al. 2015).

Armed with these landscape complexity maps, we hope to provide protected area managers with a blueprint to start disentangling the role of major ecosystem drivers. For example, are highly complex and diversifying landscapes largely driven by herbivore dynamics, disturbance events or management action? Park monitoring and research planning can be stratified using these ‘blueprints’ to begin answering these and other questions related to landscape complexity.

Cressie et al. (2009) and Lechner et al. (2012) stress the importance of accounting for uncertainty in the analysis of complex ecological data. We highlighted here the importance of accounting for spatial structure in ecological data analysis but did not assess the influence of resolution scale on analysis results. In future studies we hope to examine these results against different pixel and moving-window sizes. How these results relate to intra-annual dynamics of land surface phenology in KNP, is another interesting question for the future (Garonna et al. 2014). Additional methods that could also be explored further to investigate nonlinear species responses in ecology include Procrustes analyses or non-metric multidimensional scaling (Borcard et al. 2011) and recurrence plots-recurrence quantification analysis (Proulx et al. 2015).

Conclusion

Despite the fact that ecological components and processes in the environment have an underlying spatial structure that is locally heterogeneous, ecological regression models often employ ‘global’ techniques which assume relationships are constant over space. Using GWR models that account for spatial variation and dependencies, we were able to provide local detail on where and when physical landscape properties drive environmental heterogeneity and how this relationship changes spatially with rainfall and season. We showed that GWR is particularly valuable for ecological studies where emergent patterns are often influenced by processes interacting at different spatial as well as temporal scales (Hewitt et al. 2007).

The spatial arrangement and magnitude of model disagreement is proposed here as a measure of landscape complexity. Areas where environmental heterogeneity is not explained by stable physical landscape properties are, instead, driven by unknown complex dynamic processes. The challenge for park managers is to identify these dynamic drivers with often limited resources. Over time, maps of landscape complexity can highlight areas where physical landscape properties remain stable drivers of environmental heterogeneity or where drivers are dynamic and signify a change in the system regime. In his review, Parrott (2010) proposes landscape complexity as a key indicator of ecosystem state. Although, further research is needed on this in the context of complexity theories, our maps of landscape complexity and its trend surface may provide insight into system regime changes. Linking dynamic processes, like herbivory, fire and climate, would therefore be a logical next step to further elucidate system functioning. Until then, maps of landscape complexity can be an effective tool for targeting monitoring and research priorities to further our understanding of the drivers of environmental heterogeneity and biodiversity.