Gradient-Based Assessment of Habitat Quality for Spectral Ecosystem Monitoring

The monitoring of ecosystems alterations has become a crucial task in order to develop valuable habitats for rare and threatened species. The information extracted from hyperspectral remote sensing data enables the generation of highly spatially resolved analyses of such species’ habitats. In our study we combine information from a species ordination with hyperspectral reflectance signatures to predict occurrence probabilities for Natura 2000 habitat types and their conservation status. We examine how accurate habitat types and habitat threat, expressed by pressure indicators, can be described in an ordination space using spatial correlation functions from the geostatistic approach. We modeled habitat quality assessment parameters using floristic gradients derived by non-metric multidimensional scaling on the basis of 58 field plots. In the resulting ordination space, the variance structure of habitat types and pressure indicators could be explained by 69% up to 95% with fitted variogram models with a correlation to terrestrial mapping of >0.8. OPEN ACCESS Remote Sens. 2015, 7 2872 Models could be used to predict habitat type probability, habitat transition, and pressure indicators continuously over the whole ordination space. Finally, partial least squares regression (PLSR) was used to relate spectral information from AISA DUAL imagery to floristic pattern and related habitat quality. In general, spectral transferability is supported by strong correlation to ordination axes scores (R2 = 0.79–0.85), whereas second axis of dry heaths (R2 = 0.13) and first axis for pioneer grasslands (R2 = 0.49) are more difficult to describe.


Introduction
In response to the Convention on Biological Diversity (Rio de Janeiro, 1992), the European Union adopted the Habitats Directive for the establishment of a coherent network of protected sites for rare, threatened, or endemic species and habitat types. This network, called Natura 2000, is aimed at preserving and restoring ecological interdependencies, dispersal, and establishment processes. European Union members need to report on their conservation status every six years. It has become clear that extensive efforts are required to obtain regulatory, technical, and scientific information as well as comprehensive ecosystem management [1]. In particular, there is a need for ecological research to be carried out beyond the local scale to implement controllable management systems. To obtain relevant knowledge about the spatial dynamic of ecological processes that influence the conservation status of habitats, spatially explicit data on the location and distribution of species are required [2].
Recent developments in remote-sensing techniques have increasingly allowed for a detailed description of spatial organization of habitat characteristics and driving environmental factors [2][3][4]. However, currently, only a few studies have implemented ecological knowledge in remote-sensing-based assessment systems for Natura 2000 monitoring [5][6][7]. There is still a considerable gap in knowledge transfer between remote-sensing specialists and ecologists in conjunction with the application demands of legal authorities [5,8,9]. The first steps in combining ecological knowledge with Natura 2000 habitat management are usually carried out using indicator species mapping [10][11][12], whereby habitat types and indicator species for habitat-status assessment are modeled separately or on the basis of object classes describing habitat quality and quantity in aggregate forms as habitat units [13][14][15][16]. Such approaches start from the premise that vegetation and habitat structures exist in a discrete pattern that can be classified a priori into categories [17]. It is assumed indirectly that habitat types and conservation status can be described by co-occurring species assemblages, as stated in the concept of ecological community assembly. The basic problem of these models is that the categories depend on ad hoc hypotheses on the observed and expected ecological relevance and cannot be adapted to new findings or changes without changing the whole model. Moreover, multiple species gradients are aggregated within a limited number of categories in which derived biotope/habitat types become difficult to interpret in terms of both class membership and spectral representation [18].
There are different approaches regarding the spatial analysis of species assemblages. A number of basic concepts, e.g., distance decay and fractal scale, as summarized in Palmer and White [19], suggest the concept of vegetation continuum [20][21][22] as a more universal description of vegetation structures. It is generally stated that vegetation compositions vary continuously along environmental gradients. Fractal self-similarity of spatial vegetation pattern is solved by setting the observation scale to individual species abundances. Species assemblages are used to describe vegetation as a whole. Therein, plant species variations are capable of representing the negative relation of distance and similarity in ecological phenomena as evidence of species turnover along an environmental gradient. Transitions are no longer unexplained sources of variance. In fact, they are thought of as fundamental properties of vegetation. In particular, management strategies need to focus on these transitional ecotones, where species richness is occasionally maximized, and competition increases sensitivity on external factors [23,24]. Gradients between or at the edge of community clusters are likely to represent patterns of processes that determine habitat structure. Such multidimensional transition areas are of utmost importance in ecosystem management as required in the Natura 2000 network, where gradual differences in habitat conditions determine the required management actions [25]. In contrast to a pre-definition of discrete habitat units, n-dimensional representation of species-environmental interrelations can be described quantitatively using ordination techniques [26].
Floristic ordination spaces have been proven to be statistically coherent with spectral signatures extracted from remote-sensing images. There are several studies relating ordination-space arrangement, e.g., of heathlands [27][28][29], bogs and wet meadows [30][31][32], tree species [33], plant strategy types [34], and plant functional responses [35] to spectral gradients, whereby evidence for spatial prediction capabilities is provided. However, to date, no detailed analyses of the Natura 2000 habitat-type-specific ordination arrangement for management purposes have been published. This study was designed on an interdisciplinary basis to describe ecologically and predict spectrally the Natura 2000 habitat types and their conservation status on the basis of floristic gradients in an ordination space. We want to find out which habitat types and pressure indicators are adequately represented in ordinated structures. It is intended to reveal habitat transition as well as habitat threat owing to species shift induced by e.g., habitat management, as reflected by species gradients in a vegetation continuum. Such habitat quality parameters are required for reporting Natura 2000 conservation status in a six-year cycle. We are, furthermore, interested in determining whether habitat types and related pressure indicators can be modeled using hyperspectral reflectance signatures. Spatially explicit transfer of habitat characteristics can help to establish area-wide remote-sensing-based monitoring systems for the conservation of valuable natural habitats. Thereby, the mapping of gradual changes in plant species and habitats shall give a detailed representation of ecological interdependencies for selecting optimal management strategies. This paper introduces a methodological framework for integrating ecological knowledge into habitat conversion monitoring. It demonstrates a combined procedure of habitat conservation status assessment from a species ordination and hyperspectral image predictions. For this purpose this study is directed by three key hypotheses: (a) The floristic variety can be described by ordination; integration of new species does not change the ordination space fundamentally; (b) Habitat types, transitions, or pressure indicators can be described continuously within the specific ordination space via spatial correlation functions; on that basis a Natura 2000 habitat conservation status assessment can be derived for management purposes; (c) Distinct habitat type areas in the ordination space can be related to patterns of reflectance.
In this study, an approach is presented that reveals the transition between habitat types as well as modulations in pressure affecting the conservation status of habitats. For the first time, an evaluation of management efforts is derived directly from an ordination space, as reflected in hyperspectral imagery.

Study Area
The study was implemented on a former military training area, Döberitzer Heide, located at 53° latitude North and 13° longitude East in the west of Berlin, Germany ( Figure 1). As a result of long-term military use, open dryland assemblages established on glacial ground moraine deposits that are mainly characterized by sandy, acidic substrate in Regosol, Cambisol, and Podzol soil types (World Reference Base) [36]. Translocation of soil substrate during military actions is reflected in a small-scale floristic variability with mosaics and interpenetration of xeric sand grasslands, herb-rich grasslands, dry heath, and pioneer woods. The main area of 3946 ha is protected as a Special Area of Conservation (SAC) within the European Natura 2000 network. The SAC includes habitat types (Lebensraumtyp (LRT)) such as Inland dunes with open Corynephorus and Agrostis grasslands (LRT 2330), European dry heaths (LRT 4030), and Xeric sand calcareous grasslands (LRT 6120). Within the study area, these Natura 2000 habitat types can be characterized by major indicator species according to Zimmermann [37]. The most prevalent indicator species are Corynephorus canescens for LRT 2330, Calluna vulgaris for LRT 4030, and Festuca brevipila grouped into Festuca ovina agg. for LRT 6120. Natural succession takes place in various patterns and different phases, just as a bundle of management activities is realized in order to preserve habitat quality. Especially open pioneer stages are threatened owing to degeneration phases where cryptogams (e.g., Cladonia sp., Polytrichum piliferum) and different grass species cover increase. Within the entire area, open drylands are generally affected by scrub encroachment (e.g., Populus tremula, Sarothamnus scoparius) and the invasion of highly competitive grasses (e.g., Calamagrostis epigejos). Heathland conversion is additionally characterized by grass encroachment (e.g., Deschampsia flexuosa) and degeneration phases where mosses and lichens cover increase as the canopy of Calluna decreases [38]. Calluna heathlands are widespread over the whole study area with varying habitat quality conditions. The conservation of open pioneer stages is mostly realized in coherent areas where heathlands and different grasslands types are adjoined. The distribution of typical xeric and sand calcareous grasslands is patchier, with only rare sites reaching a good conservation status. Soil substrate variations particularly influence the quality of calcareous grassland habitats by inducing species shift along acidity gradients (e.g., Luzula campestris). Since 2004, different strategies of habitat management have been implemented by the nature foundation Sielmanns Naturlandschaften. These include the repressing of tree species or highly competitive grasses growth through big mammal grazing (e.g., Bison bonasus and Equus ferus przewalski), tree removal, and mulching of Calluna heath to support regeneration.

Floristic Data
In order to determine the vegetation continuum for open dryland habitats (including LRT 2330, LRT 4030, and LRT 6120) of the research area, vegetation samples were collected on 58 plots. Species abundances were estimated using the enhanced Braun-Blanquet method [39], whereby species nomenclature is based on Rothmaler et al. [40]. Additionally, for every plot the Natura 2000 habitat type as well as the habitat conservation status was mapped. Terrestrial mapping of conservation status was conducted using the national assessment scheme framework proposed by "Bund/Länder Arbeitsgemeinschaft Naturschutz, Landschaftspflege und Erholung" (LANA) [41] and adapted for the federal state of Brandenburg by Zimmermann [37]. It incorporates the core assessment criteria-habitat structure, species inventory, and habitat disturbance-towards three assessment categories for a favorable (A: excellent, B: good) or an unfavorable (C: adverse) conservation status. All criteria are defined by thresholds of plant species abundances and expert evaluations (e.g., present, low, extensive) [37] integrating characteristic communities of habitat conversions that are typical for our study area (see Section 2.1). Consequently, habitat pressure, represented in B/C assessment categories, can be described by structural parameter (e.g., senescence, vitality) and listed plant species assemblages [37]. Pressure strength is maximized when (a) structural and species diversity is low or (b) the influence of disturbance species is high. On the basis of expert knowledge, the spatial distribution of the sample plots was chosen so as to cover all relevant vascular plant species, mosses, and lichens, thus including all important habitats with typical transitions, succession states, and pressure indicators. In total, the fractional cover of 98 species was estimated in 1-m 2 plots. To ensure that the vegetation properties can be adequately mapped with hyperspectral imagery, the plots were located within homogeneous structures according to species composition, bare soil, and litter cover within a minimum radius of 5 m.

Species Ordination and Floristic Pattern Significance
In our first hypothesis, we argue that only a stable and significant floristic pattern, reflected in an ordination space-derived vegetation continuum, can be used to describe habitat characteristics for management purposes. We applied a nonmetric multidimensional scaling (NMS) procedure [42] on a site-by-species matrix to project rank-ordered species similarities into two-dimensional ordination space ( Figure 2A). The original number of plant species was reduced to omit species that rarely appear with low abundances over all field plots. These are known to produce strong distortion effects on the final ordination topology without increasing floristic pattern significance [43]. Furthermore, owing to a weak spatial representation, their introduced variance cannot be assumed to be causally related to image spectra. Similarities were then calculated using the Bray-Curtis distance measure [44] on the final matrix of 58 sites by 38 species. We used Kruskal's stress value [42] to interpret the goodness of fit for the resulting ordination space topology. To avoid local minima for stress values, the procedure searches within 1000 random start configurations until a stable solution is reached.
Since an ordination space for species assemblages is a generalized representation of the ecological environment, projected floristic patterns need to be assessed on their ability to represent ecological relevant structures. Furthermore, the stability of the projected patterns reveals whether an appropriate sample size was chosen to describe floristic heterogeneity adequately. Hence, we define two null hypotheses stating that there is no stable ordination plot configuration, and the ordinated pattern is not significantly different from random configurations. We used a combined statistical algorithm, testing sample stability and structural strength on ordination axes scores introduced by Pillar [45,46]. Stability was tested by generating 1000 bootstrapped samples [47,48] from the final site-by-species matrix. The bootstrapped matrices were then projected into ordination space with NDMS transformation and axes scores were compared to reference ordination after score matrix matching by Procrustes adjustment [49]. Subsequently, stability (C) was evaluated using the average Pearson product moment correlation (r) between reference scores (S) and test scores (S*) in each ordination dimension (i) over all bootstrapped samples (n): , * / .
Pattern Significance was tested, generating 1000 random permutations from the final site-by-species matrix. Permuted scores were calculated using NMS transformation and compared with test scores taken from a second NMS on the permutation matrix using the same bootstrap samples as derived in the stability test. Permutation scores (S p ) were then correlated (r) to the bootstrapped permutation scores (S**) and results were compared to the bootstrapped correlation from the stability test. We then calculated the probability (P) of permutation correlation being greater or equal to our reference correlation over all bootstrapped sample (n): We can now reject our null hypotheses for 1 − C < α and P < α, respectively, whereby α probability threshold was defined with 0.10.

Habitat Type and Habitat Pressure Aggregation
Aggregation techniques are needed in order to translate species composition of ordination plots into Natura 2000 habitat categories ( Figure 2B). On the basis of expert knowledge, site-specific vegetation characteristics (see Section 2.1), and listed Natura 2000 habitat indicator species [37], a functional plant species relation was developed for habitat type and habitat pressure evaluation. Specific habitat functions consist of a weighted sum of cover values for indicator species (Table 1). Again, weights are defined by expert knowledge incorporating site-specific habitat characteristics and legal requirements for the conservation status assessment. The weighted aggregate of habitat function components was standardized between 0 and 1 over all plots to represent a probability scale in case of habitat-type aggregates or a relative strength of influence for pressure aggregates. Standardization was performed by dividing the weighted sum of a plot by the maximum that can be reached considering probabilities in all plots. Every plot can be uniquely defined by score coordinate pairs at positions u in the ordination space. Thus, we can describe information related to plots as a realization z(u) of a spatial random variable Z that holds the distribution function for all possible realizations [50]. A realization of a habitat/pressure function can consequently be written as where β denotes the weights of the components (e.g., plant species) N for the plots i-n. Single components and related weights were selected as indicators for defining the habitat types (typical habitat indicators) as well as pressure parameters (negative pressure indicators) to assess the conservation status (Table 1). We thus assumed that the habitat indicator species would be positively linked to the occurrence probabilities of habitat types when they are known as typical character species. A negative link can be discerned when they are considered to be pressure indicators for habitat conversion. Finally, probability/strength values can be assigned to plots in the ordination space as discrete translation of the allocated species composition.

Surface Analysis and Interpolation in the Ordination Space
Our hypothesis states now that z(u) is spatially determined and therefore can be described by spatial correlation functions to predict habitat-type probabilities and pressure strength on unknown grid cells for the entire ordination space ( Figure 2C). However, as a nature of ordination, similar information is grouped in clusters with gradual changes to adjacent regions with different floristic compositions [51]. This trend violates the intrinsic hypothesis as an assumption for geostatistical prediction [52] and superimposes inner group variability that should be detected in order to assess habitat quality. To overcome this, we first separated the spatial trend. This was done by fitting first-, second-and third-order polynomial regression models for score axes with ordinary least squares (OLS). The best model according goodness of fit was selected to predict the broad scale trend of habitat type characteristics within the ordination space. Subsequently, a variogram analysis was carried out on the model residuals. We used the geostatistical approach, which combined spatial correlation modeling (variography) with subsequent spatial predictions (kriging) [53]. Herein, spatial correlation functions can be modeled by fitting an experimental variogram that describes spatial variance γ for plots i in relation to distance classes h. Every habitat function is assumed to have a typical correlation length (range) at which the maximum variance (sill) between point pairs is achieved. From that range distance, the variance decreases towards zero distance where an inexplicable minimum variance (nugget) remains. From this, one can then describe spatial correlation structures using variogram models fitting nugget, sill, and range parameters within the codomain of the spatial boundary condition of the ordination space [54]. We used an effective range in which 95% of the maximum variance was achieved to interpret the correlation lengths. Furthermore, we introduced a modified coefficient of determination, R 2 var, to describe the amount of explained variance for variogram models in comparison with a null model. As an appropriate null model where no spatial correlation could be identified, we selected the nugget effect model with no range parameter owing to maximum variance levels over all distances. The nugget level was defined as the median variance for all possible pairwise distances (sample variogram). We then built the ratio between the sum of squares for variogram model residuals (SSR) and the sum of squares for null-model residuals (SSN). According to R 2 var = 1 − SSR/SSN, spatially determined habitat functions can be identified when their variogram models contribute significantly to the explanation of spatial variance.
A list of 19 different variogram models was fitted to residuals using generalized least squares [55]. The model with the best fit regarding the minimal sum of squared error [56] for variances at all pairwise sample distances was selected to describe the spatial autocorrelation and calculate the kriging weights. Kriging was applied on a regular grid with 0.01 intervals that was expanded inside the score axes. This procedure was applied to (a) field-based habitat types and conservation status assessment and (b) habitat-function-based habitat types and pressure strength. For terrestrial habitat types, we used regression kriging of indicators [57], adding Krige interpolation and predictions from a logistic regression. A logit link function was used to transform the final results to occurrence probabilities. Simple regression kriging with a polynomial regression was applied to terrestrial habitat assessment categories and habitat-function-based habitat type probabilities and pressure strengths. In order to identify significant trend axes for regression models, we applied a backward variable selection until the Akaike Information Criterion [58] was minimized. To compare the goodness of fit for coordinate regression approaches, we used adjusted R 2 and, for a better comparison, the Nagelkerke R 2 [59] in the logistic regression models.
For external validation purposes, we compared the final variogram models and resulting kriging interpolations for both field-based and habitat-function-based derivations of habitat type and habitat pressure. To show how terrestrial mapping as reflected in ordination structures can be reproduced on the basis of functional relations that regularly connect plant species occurrences, the resulting kriging grids of both mapping methods were correlated. The average deviation between all kriging pixels was evaluated using the Pearson Product Moment correlation as well as the R 2 in a linear regression. Additionally, the variogram model parameters were compared in order to estimate the spatial correlation strength of habitat type and assessment/pressure within the ordination space.

Habitat Transition and Habitat Pressure Analysis
Isosurfaces derived from the combination of trend surface modeling and kriging predictions can be used to identify habitat type transitions and habitat pressures by means of isosurface recombination and reallocation of information stored in ordinated plots ( Figure 2D). Habitat-type probabilities are generally constructed to reveal the potential of habitat type establishment on the basis of typical habitat indicator species. To clearly demonstrate transition zones, we combined the occurrence probability grids by multiplying probabilities less than 50% for specific habitat type pairs. We assumed that below this individually replaceable threshold, ordination space can be used to reveal inter-habitat-type transition as typical habitat conversion. Above this threshold, we assumed that more distinct species-dependent pressures to habitat quality can be revealed. The strength of inter-transition is derived by a min/max normalization of the arithmetic product of probability surfaces for habitat type pairs.
Intrahabitat pressures that are responsible for the threat of habitat quality can be revealed by defining a habitat function on the basis of weighted indicator species (Table 1). Here, the relative strength of pressures allocated to a habitat type with an occurrence probability above 30% is calculated as a realization of z(u). Consequently, the strength of influence is positively correlated to the number of pressure species and their fractional cover. More specifically, areas of strong pressure influence were categorized on the basis of species compositions reallocated to distinct ordination regions. The strength of individual species influence was calculated with a min/max weighting according to specific species cover of related plot position in the ordination space.
For the purpose of conservation status assessment, we combined habitat type probability functions with pressure strength functions. We assumed that the probability of a certain habitat type is reduced when pressure factors increase. In conclusion, predicted ordination space grids for habitat type occurrence probabilities were subtracted by pressure-strength grids. The result was equally scaled to three different color intensities with gradual transitions. Finally, we categorized three assessment levels (A: excellent, B: good, C: adverse; see Section 2.2) in the center of each color class, whereas habitat probabilities ≤0% were excluded from the visualization.

Spectral Data
Hyperspectral images were acquired during a flight campaign on 4 June 2011 between 10:00 and 12:30 (Universal Time). The imaging spectrometer used was the Airborne Imaging Spectrometer for Application (AISA DUAL (UFZ, Leipzig, Germany)) ranging from visible (400 nm) to shortwave infrared (2500 nm) in 367 spectral bands. The pushbroom scanning system operated in a 24° field of view with an instantaneous field of view of 0.075° for the coverage of single ground elements. In total, 22 flight stripes with 300 samples per scanning line were recorded. The mean flight altitude was 1500 m above sea level, and the mean aircraft speed was 180 km/h. Images were geometrically corrected using an inertial measurement unit and ground control points. Overlapping flight stripes were merged into a single mosaic using an adjusted algorithm for automated control point allocation (Scale Invariant Feature Transform) [60]. The final product pixel size was resampled to 2 m × 2 m. Internal radiometric calibration was supplemented with spectral binning, smear correction, and destriping (Reduction of Miscalibration Effects) [61] to generate reliable at-sensor radiance. In order to obtain top-of-the-canopy reflectance (TOC), a radiative transfer model (ATCOR-4,) was implemented, followed by an empirical line correction (ELI) [62]. As a reference for ELI post-calibration, we used field spectra that were collected around the acquisition time with a field spectroradiometer (ASD Inc., Boulder, CO, USA). To account for observed nonlinearity within a range of 400-600 nm, we adjusted the usual ELI procedure with polynomial regression equations until the best polynomial fit between the image and the reference spectra was found. Reflectance signatures of the field plots were finally extracted from the image mosaic. A transformation to 1035 spectral variables including continuum removal [63], first Savitzky-Golay derivative [64], and spectral indices for water, pigment, nitrogen, cellulose, lignin absorption, and band-depth-normalized absorption features [65] provided spectral predictors for a coherence analysis with ordination space arrangement. The continuum was derived by fitting a convex hull over the top of a reflectance spectrum. Subsequently, absorption features are generated by dividing the original spectrum by the continuum curve. Savitzky-Golay derivatives are produced on the basis of a second-order polynomial filter of the original spectrum. The first derivative was calculated stepwise for a five-point filter length in order to render the slope for the entire spectrum. The derived spectral variables are listed in Table S2 in the Supplementary Materials (Supplementary B). Narrow spectral bands as well as overlapping physical plant properties lead to redundant spectral information. Redundancy in statistical models causes problems of multicollinearity with unreliable estimates of regression coefficients [66,67]. We therefore used partial least-squares regression (PLSR) [68], which calculates the orthogonal linear combination of original predictor dimensions (latent variables). A variable pre-selection can increase the predictive power of regression models [69,70]. Hence, dimension reduction in latent variables was incorporated with backward variable selection using a wrapper approach maximizing the model's goodness of fit based on predictor significance and variable importance implemented in the R package autopls [34]. Separate models were generated for axis scores as dependent variables. Within an internal leave-one-out (LOO) cross-validation, the number of latent variables for the best model was estimated, minimizing the error of prediction. LOO statistics were used to evaluate the predictive accuracy [root-mean-square error (RMSE)] and goodness of fit R 2 for individual axis models. Furthermore, the number of selected latent and predictor variables was used to evaluate PLSR model stability. Thereby, an increase in model complexity is incident to the consequences of model overfitting (variance-bias trade-off).

Ordination Space Stability and Pattern Significance
The final two-dimensional ordination space that showed the floristic variance distribution within our study area yielded a stress value a = 0.0016. This can be interpreted as an excellent representation of initial species composition [42,51]. The cover values of the major indicator species (see Section 2.1) are well separated into different ordination plot regions with their transitions (Figure 3a). Although a third of all samples per bootstrap iteration were excluded from the NMS ordination in each bootstrap iteration [71], the average correlation over all iterations with n = 1000 samples was high at C = 0.969 for the first axis and C = 0.956 for the second axis (Figure 3b). The interquartile range (IQR) is higher for the second axis, with more outliers to lower correlation. Nevertheless, the difference 1 − C for averaged correlations was lower than the α threshold 0.10 for both axes. Hence, we can reject the null hypothesis and state that the reference ordination space is stable in terms of plot selection. Comparing bootstrapped samples from randomly permuted data with the same bootstrap sampling units, we can see an increasing IQR with correlations ranging from 0 to 0.93 (Figure 3b). Thereby, the averaged correlation of the first ordination axis amounts to C = 0.714, and for the second axis C = 0.629. With a probability of P = 0.033 for the first axis and P = 0.021 for the second axis, the permuted correlation is higher over all iterations. Again, the α threshold was undershot, and it could be alternatively assumed that reference ordination space represents significant floristic structures.

Variography
For the three main habitat types in the open drylands of the Döberitzer Heide, we fitted variogram models to predict the occurrence probability of habitat types and the relative strength of pressure factors to assess conservation status on the basis of the habitat functions on the ordination plots. The results were compared with plot-specific field-survey data, including habitat-type delimitation and habitat conservation status assessment ( Table 2). As expected for all habitat types, a significant spatial coherence can be observed for both ordination axes, except for LRT 2330, where only the NMS2 direction features a significant trend. Comparing R 2 reg, it can be clearly seen that a spatial trend is more influential on habitat-type transition (R 2 reg Habitat type probability ≫ R 2 reg Pressure strength) for both habitat functions and terrestrial datasets, whereas change owing to pressure indicator species is more dependent on the floristic composition for LRT 2330 and LRT 6120, as reflected in higher values of R 2 vario that explain the residual variance. It can generally be revealed that species-rich plot compositions show a lower spatial dependency, which is particularly evident for LRT 6120 where R 2 reg ≪ R 2 vario. Generally, variogram models are able to explain plot variances of experimental variograms from 69% to 95% in eight of 10 cases, considering R 2 vario. Only variogram models for pressure factors and the assessment parameter for LRT 4030 are less than 50% better than a null model. In the case of LRT 2330, variogram models can explain spatial variances even better than terrestrial data. In all cases, an effective range up to a maximum variance, that is, at least 68% higher than the nugget variance, can be derived. Table 2. Variogram models for field-survey-based habitat types and habitat conservation status assessment (ter.), and for habitat-functions-based habitat types and pressure strength (fun.). Mat = Matern with kappa = 5; Cir = circular; Sph = spherical; Ste = Matern with M. Stein's parameterization; cn = nugget; c0 = sill; a0 = effective range; R 2 vario = coefficient of determination for variogram models; R 2 reg = coefficient of determination for coordinate regression; dim reg = significant dimensions (v1, v2) in spatial regression.

Habitat Type Functions and Assessment of Pressures
Using relevant habitat type functions with specific variogram models, the occurrence probability for three different habitat types was spatially predicted within the ordination space on the basis of kriging weights (Figure 4). Thereby, isolines represent locations with equal probabilities, whereby the 30% threshold of being a specific habitat type is highlighted with a dashed line in bold black. For all three habitat types, clear separations into different ordination space areas with typical inter-habitat transitions could be identified. Whereas LRT 4030 shows an omnidirectional decrease in occurrence probability, it can be shown that the distribution of habitat function components, bare soil cover for LRT 2330, and Agrostis capillaris and Festuca ovina agg. for LRT 6120, is more variable. These components overlap with adjacent habitat type distributions, whereby habitat conversion through transition is made visible. Furthermore, variations of occurrence probabilities above 50% occur as a result of varying indicator species abundances owing to the presence of pressure species. Habitat type transition within the ordination space is visualized in Figure 5. The first transition is located between pioneer stages of inland dunes and dry heath. This gradient of overlapping probabilities is mainly characterized by a change in lichen cover. The second transition between European dry heaths and Xeric sand calcareous grasslands is realized in two situations. Changing cover of different grass species on ordination plots is overlain with decreased Calluna vulgaris proportions in the upper part. A direct transition to LRT 6120 in the lower part is based on a change in characterizing herb cover. This transition is weaker because a typical herbal diversity for LRT 6120 may not be directly linked to heathland transition. The typical transition is weakened by intermediate grass stages such as Festuca ovina agg. or Nardus stricta. Within the ordination space, no direct conversion between LRT 6120 and LRT 2330 can be identified.  Figure 6 shows the kriging-predicted pressure strengths for chosen habitat types. For all habitat types, locations with strong pressure influence can be detected. In contrast, there are stable locations where there appears to be no influence of any pressure species. The plot-specific informational content within the reference ordination space can be subsequently used to assign pressure factor complexes for interpretation of habitat structures (Table 3). Regarding the habitat-quality status of LRT 2330, an important threat can be seen in a loss of bare soil cover with increased lichens and moss cover (Aa in Figure 6). This status changes into strong pressure complexes of establishing Rubus shrubs interspersed with Rumex acetosella and moss species (Ab in Figure 6). An increasing Rumex acetosella cover is also linked to an increased pressure of grass invasion (Ac in Figure 4). In particular, Agrostis capillaris cover can be identified as an important parameter for grass invasion, while its presence is often connected with xeric grassland herbs (Ad in Figure 6).
The composition of intra-habitat pressures is more complex within LRT 4030. We can discriminate between different grass invasion categories. While Ba-c in Figure 6 is dominated by a transition between Festuca ovina agg. and Calamagrostis epigejos communities, the Bd-Bf gradient in Figure 6 is characterized by Nardus stricta and Deschampsia flexuosa mixtures. These gradients are well defined at the transition to LRT 6120 and can be transferred to a better differentiation of grass invasion categories. In addition, ordination space arrangements enable the identification of shrub invasion with Sarothamnus scoparius (Bc-Bd in Figure 6) as well as tree establishment (Bg in Figure 6), which is superimposed with increased lichen cover.
Base-rich and herb-diverse LRT 6120 habitats occupy only small areas of the ordination space. These are often adjacent to grassland species that can also become established under acidic conditions. The predicted pressure strength reveals different gradients for grass species (Cc-Cf in Figure 6) that are not characteristic for a favorable status of LRT 6120. Thereby, the ordination space arrangement can be used to separate typical habitats from various different grassland types. Furthermore, pressures through tree growth in Ca-Cb will have a strong influence on habitat quality. Table 3. Pressure-complex definition on the basis of plot localization within a region of maximum pressure strengths on the ordination plane. Species cover is aggregated over a certain number of plots by min/max-normalized fractional cover values in order to assess the direction of species influence on habitat pressures.   Table 3, and dashed lines denote a habitat type probability of 30%.
On the basis of derived inter-habitat transition and intraspecific pressure complexes, a Natura 2000 habitat type assessment of conservation status was realized. This results in a grid-based continuous assessment for ordination space locations dependent on habitat type positions (Figure 7). Thereby, the conservation status can be described by three color intensities with gradual transitions. The center of each habitat type represents a favorable conservation status, whereby internal fluctuations and inter-habitat transitions are characterized by decreasing habitat qualities. We validated the distribution of conservation status assessment within the ordination space, calculating the Pearson product moment correlation and RMSE for assessment grids derived for field-survey-based assessment functions (Table 4a). Over all habitat types, a strong correlation with field surveys can be observed. The generated ordination space assessment approach differs at most by 15% from terrestrial assessment, which is within the range that can be achieved by subjective human differences. The lowest Pearson correlation with field surveys occurs for LRT 6120 (<0.859), which is also evident in habitat type prediction. In general, habitat type and assessment functions, generated by floristic composition on ordinated plot location, can adequately reproduce results obtained from terrestrial mapping in the study area. Table 4b provides a summary of the habitat-type-specific spectral PLSR model parameter and LOO accuracy assessment. Regression models that relate reflectance to scores on the first ordination axis can explain habitat-type-specific variances of up to 82% in internal validation. The lowest fit was generated at LRT 2330, where 49.1% of score variability could be explained by spectral variables. This resulted in a maximum RMSE of 21%. In second-axis models, the RMSE is maximized for LRT 4030 (RMSE = 20%). The related model provides a poor explanation for the variance in the second ordination dimension (R 2 = 0.13). In contrast, the explanatory power of second-axis models is high (R 2 > 0.80) for LRT 2330 and LRT 6120. The number of latent variables selected is small: n_C = 2 for all models. This small number of latent variables indicates model stability, owing to a high score variance, which can be explained by a minimal number of orthogonal components in PLSR. The original 1035 spectral variables were drastically reduced between 147 and 9. In particular, species-rich LRT 6120 can be explained spectrally by only a small number of significant spectral variables on the ordination plane. Table 4. (a) External validation between kriging grids on the ordination plane for terrestrial mapping and habitat functions. (b) Internal LOO validation between spectral variables and axis scores. cor = Pearson product-moment correlation; RMSE = root mean squared error; R 2 = coefficient of determination; n_C = number of latent components in final PLSR-model; n_pred = number of significant predictors/spectral variables.  In order to prove model transferability and demonstrate spatially explicit habitat type monitoring, we applied PLSR models on an open dryland area of the Döberitzer Heide. There, habitat type occurrences as well as related conservation status assessment for >30% occurrence probabilities were predicted after masking any tree and shadow pixels ( Figure 8). Generally, a clear distribution pattern of specific habitat types can be mapped. Results indicate that the typical floristic composition for habitat type LRT 6120 characterization is present in only a few regions (probability >40%). This is also reflected in predicted assessment categories where conservation status is mainly assigned between C and B (unfavorable).

Spectral Predictability
In fact, habitat type LRT 6120 occurs in various transitions to pioneer grasslands and dry heaths as shown in red (10-40% probability). Open pioneer grasslands and dry heaths are more common in the study area. Their conservation status mainly ranges between A and B, whereas spatial patterns indicate an expected decrease in habitat quality from core areas to edge regions (Figure 8 zoomed subplots). External validation was performed on the 58 field plots by extracting habitat types for a probability threshold of >30% and for equally spaced assessment categories. Habitat types LRT 2330 and LRT 4030 can be mapped with an overall accuracy (OAA) of 100%, whereas species diversity in LRT 6120 is more difficult to detect (OAA = 73.3%). However, degenerated stages of LRT 6120 with probabilities <30% were not included in the validation. Terrestrial assessment categories show good conformity with LRT 2330 (OAA = 84.2%) and with LRT 4030 (OAA = 89.5%). Conservation status variations are more complex in LRT 6120, which results in an OAA of 66.6%. ; a typical transitional area between the three habitat types was exposed in the subplot zoom.

Spatial Correlation
Our study demonstrates the use of spatial correlation functions to determine habitat types, pressures, and conservation status in a site-specific ordination space. As an initial step, we introduced habitat functions as representations of habitat occurrence and pressure/threat strength. It should be noted that predicted habitat patterns are strongly dependent on selected species and chosen species weights. In this respect, our study presents a straightforward procedure to determine how expert knowledge on habitats and habitat pressures can be transferred to ordination space projections. The modeled type and status therein are seen as possible representations of ecological interdependencies in a vegetation continuum. There is no general allocation of floristic composition to a certain habitat type or pressure complex. Every ordination space can be quantified individually according to the study area, assessment demands, or management purposes. Our approach provides a reproducible aggregation technique on the basis of species lists and is therefore distinct from a priori habitat classification or obviously subject-dependent terrestrial assessment.
The species composition used in this study to describe the conservation status categories for dry heath is based on the legal standards defined in Annex I of the European Habitat Directive [72], as well as expert knowledge [37,73]. However, the proposed methodology is not restricted to Natura 2000 habitat types. With an appropriate sampling of indicator species and pressure factors, every monitoring or assessment approach can be analyzed on its ability to reflect clear patterns in an ecological gradient space. Thereby, habitat type probabilities as well as pressure strength are spatially predicted on the basis of variogram models. In geostatistics, there is no standard methodology to select an appropriate model. In our study, the best model was selected by minimizing the prediction error for a choice of 19 known models. Nevertheless, it is important to keep in mind that the final results for a grid-based probability pattern are dependent on the choice of spatial correlation function (variogram model) and its overall predictive capacity. Spatial probability patterns are therefore not deterministic and can only be approximated, taking into account adapted selection algorithms [74,75]. Another source of spatial uncertainty is in the ordinary kriging procedure itself. The number of points used to calculate weights for an unknown grid cell can have an influence on spatial heterogeneity. We constantly used half the number of total plots per grid cell to derive reliable kriging weights.

Species Composition
Probability aggregation in ordination space dimensions is usually applied on external variables to interpret abstract gradients. Vegetation ecologists are well aware of spatial statistic methodology [76], which is used to produce isolines representing external correlation structures by means of classification approaches [77] or trend-surface analyses [78]. To our knowledge, this is the first time that multi-species probability estimation, on the basis of habitat/pressure functions, has been examined. In addition to habitat type and threat, the conservation status can consequently be described by ordination space structures that reveal species gradients on the basis of pressure definition. However, even though separation of general gradient patterns with axis models reveals fine-scale floristic heterogeneity that can be described by means of variography, the identification of unique species complexes may become complicated in species-rich continua. Therein, differential species contribute at different gradient positions to habitat quality and distribution. Species complexes are not generally separated in single positions in ordination owing to overlay and indifferences as part of the unexplained variance. Even if habitat types can be directly allocated using probability thresholds, a distinct separation of near-ordinated but floristic variable plot locations should be reviewed critically.
Besides gradually changing species cover in adjacent plots, abrupt changes in species representation as revealed by pressure complexes (Table 3) are evident in ordinated species composition. In addition to axis stability and pattern significance estimation, a good floristic representation can be further increased by optimizing preserved sample variance. In our study we used a two-dimensional NMS with an excellent representation of the floristic variation (stress = 0.0016) in order to demonstrate a two-dimensional Kriging procedure. The decision was based on evaluating the strength of spectral correlation to single score axes. The averaged R 2 of the first NMS axis over all habitat types was maximized in a 2D solution. However, additional variance patterns may be related to spectral signatures. For this purpose, a case-specific choice of number of ordination axes, distance metric, original dimensions (surface and vegetation structure parameters besides plant species), and a detailed analysis on recent algorithmic developments such as Isomap [28,79] still ought to be considered.

Spectral Application
The spectral discrimination of axis gradients varies for specific habitat types and selected axes. It should be noted that as part of the applied NMS ordination, axes are principal component rotated in order to explain the maximum variances in the plot configuration. The resulting directions are not automatically related to spectral diversity, and it can be assumed that linking the spectral discriminability to axis-specific rotation angles will increase the predictive accuracy. Further research is needed to find supporting evidence for this. Another source of unexplained regression variance can be seen in the representation of the spectral sample itself. Spatial heterogeneity on 2 m pixel size can introduce an increased signal variance owing to adjacent effects. Furthermore, spatial non-stationarity due to phenology shift or varying litter cover can influence model representation on image pixels [16]. In addition, image-spectra calibration always delivers spectral response models under the boundary conditions of acquisition time. Spectral library information on the basis of TOC reflectance can be considered to be an improvement for true variance estimation and transferability when phenological phases are covered adequately. Nevertheless, the transferability of regression models for floristic patterns still remains complicated owing to vegetation status, irrespective of the species [80]. Additional parameters such as the chemical constituents under the influence of plant stress and growth [81], and spatial heterogeneity such as litter content and canopy height [82], should be described in order to obtain reliable models for monitoring purposes. However, the approach presented here can enhance a Natura 2000 habitat assessment with spatially explicit predictions of conservation status incorporating floristic compositions along ecological gradients.

Conservation Status Assessment
The presented approach demonstrates a pixel-wise conservation status assessment on the basis of Natura 2000 habitat type transition and pressure indicators that are directly derived from ordination space structures. An important advantage can be seen in the decoupling of the spectral and the ecological models. We can spectrally predict the vegetation continuum and a posteriori derive information from that. It crucially differs from common remote sensing based methods, where image pixels are classified according to different habitat types [12,83] or habitat quality parameters [12,15]. Therein, an image pixel is determined by one attribute that was a priori defined as a relevant ecological entity for the evaluation of habitat quality. Various habitat quality indicators are developed [10] that allow a fine-scale prioritization of management strategies. Although remote-sensing-derived habitat quality maps show a good correlation to terrestrial mapping approaches, they can only explain variations in fine-scale conservation status indicators up to 39% [6]. In the proposed approach, the information mapped at the pixel scale is variable. Fine-scale variations are directly transferred from ordination space via spectral coherences. Both habitat transition and pressure species complexes are transferable to images using the informational content of the ordination space that projects the floristic variation in an environmental space. This enables additional conclusions about the mapped conservation status. Thereby, an image pixel is linked to the structure of the site-specific ordination space that holds information such as the direction of habitat succession, the distribution of plant species, or, indirectly, about the abiotic gradients. Commonly, this information has to be defined before mapping and the conservation status assessment is based exactly on these defined categories [12,14]. The functional aggregation technique coupled with probability, pressure strength, and assessment predictions also allow a continuous interpolation of ordinated plot information. Hence, habitat conversion can be made visible in continuous gradients when ordination dimensions are transferred to image data. Development tendencies with regard to species shifts can be revealed in these transitional areas. However, the aim of the study was not to give a complete conservation status assessment. It is rather aimed at providing a methodological framework for the evaluation of plant species shift that is assumed to be responsive to management in our study area (grazing, mulching, species removal). We do not include additional, structural vegetation parameters (e.g., vitality, senescence) or anthropogenic influences (e.g., burning, nutrient transfer) in the ordination that can increase the accuracy of habitat quality assessment. It has further to be mentioned that this study is based on a site-specific ordination space for open dryland habitats on former military training areas in Brandenburg. In order to reveal fine-scale variations in transition and species composition, such ordination results are restricted to certain biogeographical regions. Integration of different habitat types always depends on the availability of species data whereby comprehensive data archives such as spectral libraries can be used to transfer the proposed methodology. Plant species data as well as related spectroradiometer measurements used in this study were therefore stored in a freely accessible database called SPECTATION [84]. Therein, field plot-specific plant species lists, vegetation class and conservation status units, and surrounding soil properties are stored for open dryland and wetland habitats in conjunction with spectral reflectance signatures for the years 2008 to 2011. This enables reproducible research on similar habitats or methodological extensions to different habitat types, which could be a subject for analysis in future studies. Species ordination and subsequent spectral variance estimation in a broader scale (e.g., country-or Europe-wide) has still to be investigated by means of new multidimensional interpolation methods. With regard to this, the crucial question for further research is: how many habitat types can we integrate in one ordination in such a way that fine-scale variations are still visible in ordination as well as in the spectral response? New statistical approaches from big data analysis in conjunction with spectral library information open future perspectives on detailed Natura 2000 habitat mapping.

Conclusions
The probability of a habitat being of a specific type depends on the habitat status incorporating inter-habitat transitions and pressure factors. The information content of ordination spaces can be used to continuously determine such habitat structure parameters. It can be shown that floristic patterns projected in the ordination space are significant and stable. There is strong evidence that functionally aggregated habitat characteristics on the basis of plant species data are spatially determined over distinct regions of the ordination space. Empirical score axes models as well as residual variogram models can be used to describe the ordination space variability of habitat characteristics such as habitat type and habitat pressure. A subsequent model combination further allows a spatially continuous interpolation of habitats and related pressure strength over the entire ordination space. Habitat transition as well as pressure indicators can be made visible in distinct ordination space regions for conservation status assessment. Results correspond well to terrestrial Natura 2000 conservation status assessment. Using evidence on spectral coherence, habitat status probabilities can be used directly to produce spatially explicit maps. This approach differs crucially from conventional remote-sensing-based habitat assessment methods that assume discrete management units as predefined natural components. Spatial monitoring is no longer dependent on threshold-based changes in habitat categories. The potential of change can be directly projected over probabilities in ordination spaces, and assessment tendencies are directly transferable to spatial information. This enables the Natura 2000 monitoring to assess habitat type vulnerability more rapidly and allows a more effective prioritization of management activities to preserve a certain conservation status. This is especially true in open land habitats on former military training areas, where habitat conversion is driven along successional gradients and terrestrial mapping is complicated by undiscovered military munition.