Organic Matter Modeling at the Landscape Scale Based on Multitemporal Soil Pattern Analysis Using RapidEye Data

This study proposes the development of a landscape-scale multitemporal soil pattern analysis (MSPA) method for organic matter (OM) estimation using RapidEye time series data analysis and GIS spatial data modeling, which is based on the methodology of Blasch et al. The results demonstrate (i) the potential of MSPA to predict OM for single fields and field composites with varying geomorphological, topographical, and pedological backgrounds and (ii) the method conversion of MSPA from the field scale to the multifield landscape scale. For single fields, as well as for field composites, significant correlations between OM and the soil pattern detecting first standardized principal components were found. Thus, high-quality functional OM soil maps could be produced after excluding temporal effects by applying modified MSPA analysis steps. A regional OM prediction model was developed using four representative calibration test sites. The MSPA-method conversion was realized applying the transformation parameters of the soilpattern detection algorithm used at the four calibration test sites and the developed regional prediction model to a multi-field, multitemporal, bare soil image mosaic of all agrarian fields of the Demmin study area in Northeast Germany. Results modeled at the landscape scale were validated at an independent test site with a resulting prediction error of 1.4 OM% for the main OM value range of the Demmin study area. OPEN ACCESS Remote Sens. 2015, 7 11126


Introduction
Owing to the demand for qualitative and quantitative soil information at multiple scales in precision agriculture (PA) [1][2][3][4][5][6], accurate functional soil maps are needed for the concept of site-specific management (SSM) to balance a profitable and cost-effective crop production with environmental concerns and sustainability [7][8][9][10].To evaluate site-specific spatiotemporal variable soil properties at the field-scale, many functional soil mapping models and approaches based on digital elevation models (DEM), proximal and/or remote sensing (RS) data have been designed.Compared with geostatistical and spatial interpolation methods (e.g., Kriging procedures, fuzzy clustering algorithms) based on comprehensive, time-consuming, and costly field surveys and soil sampling, non-invasive remote and proximal sensors, combined with empirical-or physical-based data analysis, offer potentially more effective, quick, and cost-efficient continuous direct or indirect data on physiochemical soil characteristics, which are determined by spatial, temporal or spectral sensor resolutions [1,6,[11][12][13][14][15][16][17][18].Due to the advantages of existing, easily accessible data archives, relatively low-cost and high-temporal, high-spatial resolution multispectral imagery and time series are available for qualitative and partly-quantitative soil information extraction, deduction of soil patterns, and mapping of SSM zones and soil surface units [7,15].
For the latest state-of-the-art science, Ge et al. (2011) [7], Mulla (2013) [8], Mulder et al. (2011) [15], and Panda et al. (2010) [19] provide research summaries of remote sensing systems and applications in PA and soil science (e.g., crop yield, biomass, soil organic content, and soil texture).Numerous studies have investigated the relation between soil reflectance characteristics (e.g., soil color and soil brightness) and soil properties, particularly organic matter (OM), soil moisture, and texture [1,[20][21][22][23][24].As proximal and remote sensing methods, diffuse reflectance spectroscopic techniques in the visible, near-and shortwave infra-red range were applied for OM measurements [25,26].During the attempt to determine the optimal content of soil organic matter of agriculturally-used soils, Wessolek et al. (2008) [27] pointed out that optimal organic matter supply is generally site-specific and strongly influenced by climate (e.g., temperature, precipitation), soil texture and land management/land use.Still, no OM simulation model exists which is applicable to all cultivation systems and agrarian used soils in Germany.
A monotemporal analysis with a single remote sensing dataset always contains heterogeneous spatial reflectance patterns due to the different land uses and vegetation conditions without temporal dimension.Regarding agricultural land use, static soil patterns persist, aside from temporal patterns characterised by crop types, vegetation phenology and land management practices.To reduce the disturbance of the temporal patterns, a multitemporal approach with the advantage of potentially higher pattern stability is needed.Blasch et al. (2015) [1] proposed an analysis workflow for OM estimation at the field-scale based on multitemporal soil pattern analysis (MSPA) using a multispectral RapidEye satellite time series, which was tested on a single demonstration field of the TERENO North-Eastern German Lowland Observatory (TERENO-NE) in Demmin, Northeast Germany.For the extraction of spatially-precise static soil patterns, it was shown that multitemporal analysis (R 2 0.81) of RapidEye data is more suitable than monotemporal analysis (R 2 0.73-0.79),whereby the latter strongly depends on the optimal time and condition of image acquisition.As soil-pattern detection algorithm standardized principal components analysis (PCA) was applied to preselected bare soil images.Statistical regression analysis showed strong correlation between first standardized principal component (PCst1) and laboratory-analyzed OM values.The underlying principle between PCst1-detected soil patterns and OM is that due to image pre-selection, PCst1 stores the reflection properties of bare soils in terms of albedo.Accordingly, an increase in OM content is correlated with a decrease in the albedo, with respect to spectral reflectance values of PCst1 with the consequence that the darker the soil pattern (lower PCst1 values), the higher the OM content, and the lighter the soil pattern (higher PCst1 values), the lower the OM content [1].However, high variance between low OM values and PCst1 was observed.Soil spectral reflectance could be additionally influenced by other soil properties (e.g., CaCO3, clay), even when no significant relationship between PCst1 and the parameters existed.The multitemporal characteristic allowed for the evaluation of the spatiotemporal variability in reflectance patterns.Using OM soil sampling data with the identified static soil patterns instead of all existing (temporal and static) soil patterns, improves the initial OM prediction model and, consequently, the OM functional soil map.
In this research the objective is to develop a multi-field landscape-scale version of Blasch's et al. (2015) [1] field-scale based MSPA model for OM (percentage weight loss on ignition) estimation based just on optical satellite data reflectance (here: RapidEye time series) as well as RS and GIS spatial data modeling to secure the idea of method transferability to other study areas at low costs.Therefore, it is necessary to prove whether OM explains the soil reflectance pattern.The steps to reach this research objective are to (1) evaluate the quality and robustness of the field-based MSPA to predict OM for single fields and field composites with diverse geomorphological, topographical, and pedological characteristics of the young morainic soil landscape; (2) create a representative regional OM prediction model; and (3) convert the MSPA method from the field scale to the multi-field landscape scale using the regional prediction model.The capability of the proposed MSPA-method at the landscape scale is investigated using four calibration and one validation test sites by applying the transformation parameters of the soil pattern detection technique and the best regional prediction model to a multitemporal bare soil image mosaic of all agrarian fields of the study area.For model calibration and validation 731 soil samples were collected at five test sites (composed of 16 fields) in the Demmin region, Northeast Germany.

Study Area
The 834-km 2 study area (upper left corner: 54°3′N, 12°53′E, lower right corner: 53°48′N, 13°20′E) is located near the city of Demmin (53°54′N, 13°3′E) in the northeastern lowlands of Germany in Mecklenburg-Western Pomerania (Figure 1).During the last Pleistocene period (Weichselian Glaciation), this region was formed by recurring glacial processes.Characterized by typical glacial morphological features, such as vast ground moraine areas, push (terminal) moraines, lake basins, glacial valleys, eskers, kettle holes, and sandars, the study area is representative of the young morainic soil landscape of Northern Germany [28].Furthermore, lakes (e.g., Kummerower Lake) and river systems (e.g., Peene River, Tollense River, and Trebel River) are tightly connected with the near-surface ground water table.With an altitude range between 0 m and 96 m above sea level (mean 19 m, standard deviation 14 m) and a slope angle between 0° and 61° (mean 1.1°, standard deviation 1.5) [29], the relief is relatively flat in the north and undulating/hilly with steep slopes in the south.Due to the significant changes in the relief, the parent substrate material, and the distance from the groundwater table, the soil types are spatially variable.The flat and marginal undulating plain is composed of sand-rich areas and glacial-till areas.Cambisols, Luvisols, and Albeluvisols are predominant in the sand-rich area.On the less sandy but more clayey and loamy glacial till, Luvisols, Albeluvisols, and Stagnosols evolved.In hilly terrain, these soils are frequently truncated, and colluvial soils developed from eroded loamy material or deposits over former bogs [28,30,31].Peaty soils are found in the floodplains.According to the Teterow weather station (34 km west of Demmin), the long-term (1981-2010) mean temperature is 8.7 °C and the mean precipitation is 584 mm/year; these records characterize the present climate of the study area [32,33].This area is part of the sparsely populated (69 inhabitants per square kilometer), intensively used agricultural state of Mecklenburg-Western Pomerania, which consists of 80.3% farmland, 19.5% grassland/pastureland and 0.2% other agrarian land use.In 2012, cereal crops (55.5%), feed crops (19.3%), oleaginous fruits (18.6%), root crops (3.6%), legume crops (0.4%), and others (2.6%) were cultivated [34].In general, the study area is located within the TERENO-NE, which is managed by the German Research Centre for Geosciences Potsdam (GFZ) [1].For the method development, calibration, and validation, five test sites that consisted of 16 agrarian fields (size: 7.53 km 2 ) were installed to represent the topographical, pedological and geomorphological characteristics of the study area (Figure 1).The test sites of AOI1 Jarmen (~330 ha, six fields), AOI2 Sassen (~124 ha, one field), AOI3 Buchholz (~162 ha, seven fields), and AOI5 Borrentin (~107 ha, one field) belong to the soil landscape sub-unit of "soils of ground morainic deposits and (mostly) loamy terminal moraines germination".AOI4 Beestland (~29 ha, one field) belongs to the sub-unit of "soils of sandars, dry valley sands, sandy plates and sandy terminal moraines" [28].Due to their specific topographical and geomorphological features, the AOI1 test site can be characterized as flat ground moraine (till) with a large (clayey) low depression and few kettle holes.AOI2 is characterized as undulating ground moraine (till, sand over till, and sand) with eroded hilltops, few low moor depressions and eskers.AOI3 is characterized as undulating curved push (terminal) moraine (till, locally clayey marl) with few kettle holes and partly eroded very steep slopes up to 11.5°.AOI4 is characterized as a transitional area from very flat ground moraine (till) to a former glacial valley (sand over till) with nearby drainage ditches and the Trebel River.AOI5 is characterized as undulating ground moraine (till) with one large flat depression, few kettle holes and low moor areas (partly covered by colluvium).In contrast to AOI1 and AOI3, very homogenous soil texture characteristics are mapped for AOI2, AOI4 and AOI5 (Figures 2 and 3).Owing to varying demand of remote sensing data handling provoked by different crop types, land uses, and/or strong relief influence, the test sites had to be categorized into (i) single fields and (ii) field composites.Due to consistent field management and crop type, AOI2, AOI4, and AOI5 are classified as "single field" (Figure 2).Although AOI1 and AOI3 are composed of single fields, these test sites are considered primarily as "field composites" due to different crop types and field management (AOI1, AOI3) as well as strong relief gradient (AOI3) (Figure 3).

Baseline Data
Since the aim of the MSPA workflow is to be transferable to other study areas at low costs, it is based just on optical satellite data reflectance.Thus, DEM data were not used because high-resolution DEM data are very costly and not available everywhere.Further, freely available DEM data (e.g., SRTM) would not significantly improve the modeling because of their low resolutions.In addition to using multispectral satellite imagery as model input data, extensive soil sampling and analysis were undertaken for model generation and validation, as well as data interpretation.Table 1 demonstrates the specific characteristics of all baseline data.During the dry conditions of July and October 2013, comprehensive soil sampling that consisted of 731 soil samples was performed at all 16 test site fields.Using the information on existing soil maps, DEMs and visible soil color patterns in the applied grid sampling design, a representative soil sample dataset was obtained.All sample positions were measured using a GPS with 5 m accuracy.This inaccuracy was reduced by generating a 5 m buffer around each sampling point, representing the position-specific mean values of the digital data (e.g., PCst1).For posterior statistical analysis, these mean values were extracted.Further characteristics of the sampling design are shown in Table 2.A representative soil surface mixed sample (500-700 g) composed of five subsamples within a 1 m radius and a depth of 10 cm was obtained at each location.The soil physical laboratory of the Department of Ecology of the University of Technology Berlin conducted the physical and chemical soil sample analysis: (1) particle size distribution using a combined sieve (fraction 0.063-2 mm) and pipette (fraction < 0.063 mm) method [35]; (2) soil acidity (pH) with 0.01 m CaCl2 [36,37]; (3) OM content using the loss of ignition method [38][39][40]; and (4) (calcium) carbonate content (CaCO3) [41].The particle size distribution results were applied to determine the soil texture classes according to the official German Soil Survey Description [42].

Remote Sensing Data
For the field-based transfer and for the further development of the landscape-scale MSPA [1], the same multispectral, temporally-and spatially high-resolution, state-of-the-art sensor RapidEye was used.Overall, 54 radiometrically-calibrated and georeferenced RapidEye scenes (Level 1B, Level 3A) from April 2009 to May 2014 were acquired from the RapidEye Science Archive (RESA) (Table 1).For this period of time, OM degradation effects are assumed to be negligible for the study area, because neither heavy precipitation (namely, erosion) events nor profound land use changes (except the conventional crop rotation) occurred.
All RapidEye datasets were pre-processed, i.e., atmospherically-corrected using ATCOR for ERDAS IMAGINE, resampled to a pixel size of 6.5 m, image-to-image co-registered using an in-house algorithm [43], and clipped to various extents (e.g., single field, field composite).The resulting RapidEye subsets, correspondingly-derived NDVI images (including descriptive statistics) became the basis for further image analysis.

MSPA Workflow
As this research is built upon the methodology of Blasch's et al. (2015) [1] field-scale based MSPA workflow, a brief description of the analyzing steps is given (Figure 4).(1) Selection of bare soil images in the RapidEye time series (data mining) using an automated classification based on NDVI thresholds (NDVI-MEAN, NDVI-SD) and phenology data to avoid effects of temporal pattern, namely, to identify datasets with static soil pattern; (2) Soil reflectance pattern detection using standardized PCA on all bands of selected bare soil images (image mining) to re-arrange information according to the assumption that the percentage area of change is relatively smaller than the percentage area of no change, stable areas will appear in the lower-order PCs, and changing areas in the higher-order PCs [44,45].This implies that the reflectance information of stable bare soil areas will show up in the lower-order PCs due to the pre-selection of bare soil images; (3) Evaluation of spatiotemporal soil pattern stability using statistical per-pixel analysis (stability analysis) to create an overall stability image calculated by the arithmetic mean (MEAN) of the pixel-wise reflectance change (namely, standard deviation, SD) between all identified bare soil images.A lower MEAN-SD per pixel corresponds to a more static reflectance pattern; (4) Functional soil mapping based on static soil pattern after stepwise exclusion of temporal effects applying percentile analysis to the overall stability image.For functional OM soil map generation the regression equation of PCst1 and OM values within static soil pattern is applied to the PCst1-image using the Raster Calculator of ArcGIS Spatial Analyst.

Application of Field-Based MSPA
The analysis steps of the field-scale MSPA [1] (Chapter 2.3) were applied to each field, each test site and the calibration test site composite "Demmin" (AOI_CAL: Composed of test sites AOI1, AOI3, AOI4, AOI5; namely, 15 fields) with the following modifications.
The automated classification based on NDVI thresholds (NDVI-MEAN, NDVI-SD) and the phenological situation per crop on the acquisition date [1] were modified, because of the unavailability of vegetation (field/crop) data for all test sites, to "just" using these NDVI thresholds to separate the static soil pattern and the temporal vegetation/land management pattern.After removing all cloudy/shady subsets from the time series, the resulting RapidEye subsets were automatically classified into class 1 images with bare soils/ripe vegetation (NDVI-MEAN ≤ 0.2, NDVI-SD ≤ 0.07) and class 2 images with vegetation/management effects (NDVI-MEAN > 0.2, NDVI-SD > 0.07).By iterative visual image analysis with image-enhancement techniques (e.g., contrast stretching and false-color composites), the most suitable (bare soil) images in class 1 were identified.
Furthermore, for each field, each test site and the calibration test site composite, all of the identified bare soil subsets were stacked to one corresponding multitemporal bare soil image.To detect the soil reflectance pattern, a standardized PCA based on the field-specific and test-site-specific correlation matrix was performed on each multitemporal bare soil image.Next, to analyze the average spatiotemporal stability of the bare soil reflectance pattern, the overall stability image (MEAN-SD) for each multitemporal stacked image was calculated using statistical per-pixel analysis [1].
The existence of a non-linear relationship between OM and PCst1 of bare soil images was shown at the demonstration field "Borrentin" [1].To statistically analyze this dependency, the position-specific mean values of (1) PCst1 of every multitemporal bare soil stacked image; (2) the overall stability; and (3) each soil property data (Table 3) were extracted at each soil sampling coordinate.Subsequently, after log10-transformation of the PCst1 and OM values, the non-linear relationships were computed as simple linear regression using the R-script "Fitting Linear Models".Moreover, the coefficient of determination (R 2 ) was calculated as prediction quality criteria as well as the root mean square error (RMSE) to express the model robustness and the efficiency of the prediction method by the leave-one-out cross-validation [46].A lower RMSE value corresponds to a more robust prediction model.
To evaluate the outlier-influence on the prediction model, a visual model diagnosis using the scale-location-plot were performed.Detected potential outliers were visually verified with corresponding spectral and soil property data.Owing to unclear sample assignment at the stage of laboratory and/or field work, in total 29 true outliers (4% of all samples) were identified and excluded from further prediction modeling.Based on the overall stability image, the iterative, percentile-based statistical analysis from the MSPA workflow was performed to define the best threshold between the static and temporal reflectance patterns.All OM sample data below the best-separating threshold were applied to the prediction modeling as the values that represent the static soil reflectance pattern.For functional soil mapping purposes of the OM content, the resulting regression equations were adapted to the analogous multitemporal PCst1-image.The very slightly acidic topsoil layer (pH-mean: 6.6) of all test sites has a low amount of CaCO3 data (10.9% of all soil samples) with low CaCO3 values (Mean 2.5% and SD 2.5%) (Table 3).Furthermore, 98.9% of the analyzed soil samples belong to the following main soil texture groups [42]: (1) sands (54.6%) composed principally of loamy sands (43.8%) and silty sands (10.7%) and ( 2) loams (44.3%) composed of sandy loams (28%), normal loams (15%), and clayey loams (1.4%).In total, loamy-sandy topsoil layers (71.8%) were predominant; 69.1% of all soil samples have moderate OM content (2-4 OM-%), followed by 21.2% with strong OM content (4-8 OM-%), 5.3% with low OM content (1-2 OM-%), 2.3% with very strong OM content (8-15 OM-%), 1.2% with extreme OM content (peaty mineral soil) (15-30 OM-%), 0.5% with peat (>30 OM-%), and 0.3% with very low OM content (0-1 OM-%).The main OM value range of the study area Demmin is from 0 OM-% to 15 OM-% (98.2% of all soil samples).

Development of the Multi-Field Landscape-Scale Version of MSPA
For model transfer to the landscape scale of the study area, the calibration test site composite "Demmin" (AOI_CAL; 525 soil samples) and its previously calculated MSPA results (Section 2.4) were used as a representative basic dataset for the young morainic soil landscape of the Demmin region covering the diversity of the geomorphological, topographical, and soil characteristics.To secure the identical PC-transformation parameters and consequently to enable the transfer of identical soil pattern detection characteristics from the field scale to the landscape scale, the correlation matrix of AOI_CAL was extracted at the PC-transformation step.After exclusion of all previously-identified outliers the log-transformed linear regression model was computed and consequently used as a regional prediction model.AOI2 with 147 soil samples (~22% of all samples without outliers) represents the primary validation test site.
To avoid the persistent problem of heterogeneous spatial reflectance patterns in single remote sensing datasets, for model transfer to the landscape-scale the multitemporal approach was obligatory to create bare soil image mosaics of the study area, which can be used as MSPA input datasets to cover all agrarian fields of the study area without temporal patterns.Heupel 2013 [47] developed a selection algorithm to identify the most suitable RapidEye dataset for soil surface analysis within the RapidEye time series for each field segment of a polygon shapefile.Threshold values were used for excluding negative effects, such as clouds, shadows, and vegetation (NDVI-MEAN, NDVI-CV).The product was a multitemporal bare soil image mosaic for each polygon/field segment at various times of the acquisition.
The vegetation thresholds of the in-house developed algorithm were modified using the same threshold values of the MSPA workflow (NDVI-MEAN, NDVI-SD; Section 2.4).Due to the quantity of field polygons, visually analyzing the negative effects of ripe vegetation in each polygon would be very time consuming, inefficient, and inaccurate.Here, further development of the in-house algorithm regarding the separation between ripe vegetation and bare soil might be required.Nevertheless, as the best existing option (to our knowledge), this algorithm was applied to compute bare soil image mosaics for the study area based on the shapefile of all study area field segments.
To maintain the multitemporal approach, it was necessary to compute three mosaics of the study area and substitute the three bare soil images of AOI_CAL in the field-based MSPA workflow because one mosaic substitutes for only one bare soil image.To avoid repetitions of already-used bare soil images per polygon in a previous composite, the identified bare soil images per polygon were sorted from low NDVI-SD to high NDVI-SD.The first mosaic was all bare soil images of the lowest NDVI-SD, the second mosaic was the second lowest NDVI-SD images and so forth.Therefore, three bare soil image mosaics were computed based on the three lowest NDVI-SD values of each field segment without repetition of previous NDVI-SD values.The described step to generate multitemporal bare soil image mosaics replaces the data-mining step of the modified MSPA workflow (Section 2.4).All following analysis steps of the modified MSPA workflow were applied to the three bare soil image mosaics of the study area, and functional OM maps based on the improved regional prediction model were generated.At the image mining step using standardized PCA, the previously extracted correlation matrix was used to transfer the identical transformation parameters from AOI_CAL to the study area bare soil image mosaics.

Relationships between Soil Properties and PCst1 of Multitemporal Bare Soil Images
To understand the general potential and the method transferability of MSPA, for each single field (16 field polygons), each side-by-side field composite (two composites), and one distant field composite (AOI_CAL), the relationship between the PCA-detected soil reflectance patterns (values of PCst1) and the laboratory-analyzed soil parameters (OM, clay, silt, and sand) were statistically investigated using regression analysis (Tables 4a and 5a).Exploiting the advantage of multitemporal data, the OM prediction model based on static soil pattern was calculated after the elimination of temporal influences (Tables 4b and 5b).
In general, the potential of MSPA for OM estimation is revealed for both single fields and field composites independent of the sample size/density, the number of bare soil images, and/or physical-geographical and land use characteristics.

Single Fields
Table 4a "Evaluation of MSPA" shows the number of detected bare soil images, the number of identified soil sampling outliers, the number of remaining (used) soil samples after outlier elimination, and the coefficient of determination (R 2 ) for each regression analysis between PCst1 and the corresponding soil parameter (OM, clay, silt, and sand).Further, the characteristics of the improved OM prediction model based on MSPA, such as the best percentile threshold for eliminating temporal effects, the number of sampling points used for the prediction model generation, model-quality and model-robustness criteria (R 2 and RMSE), and the results of the OM prediction model, expressed as the extremes of the modeled OM values are demonstrated in Table 4b.For all remaining OM samples, significant (p-value < 2.2 × 10 −16 ) linear (after log-transformation) regression models for PCst1 and OM with a coefficient of determination R 2 ≥ 0.66 was obtained for each single field."An increase in OM content is correlated with a decrease in the spectral reflectance value of PCst1, and consequently, the darker the soil pattern, the higher the OM content, and the lighter the soil pattern, the lower the OM content" [1].Spielvogel et al. (2004) [26] and Konen et al. (2003) [48] found a similar non-linear correlation between OM and soil lightness resp.soil color (Munsell value, chroma, and reflectance).The only exceptions are AOI3d and AOI3g: Due to relief-shadow effects (AOI3d) and a very low sample size of nine (AOI3g), no correlation was found.The problem of relief-induced shadow effects is very common in remote sensing.Since shadow areas always contain disturbed reflectance values and lead to incorrect data interpretation, they should be masked and excluded from further modeling.Subsequently, these single fields were excluded from the field-composites (AOI3, AOI_CAL) and their spatial/statistical analyses.
At the single field polygon level, the coefficients of determination (OM, PCst1) vary substantially (0.66 ≤ R 2 ≤ 0.91), with the lowest for AOI4 and the highest for AOI1c.The R 2 of all single fields of the field-composites ranges between 0.74 and 0.91 for test site AOI1 and between 0.68 and 0.81 for AOI3 (Table 4a).Note that low sample sizes and/or input data quality, besides physical-geographical characteristics (e.g., relief), can lead to fluctuations.Therefore, excluding low-sampled (<25 samples) single fields from further discussion, for some single fields, the clay (AOI1b: R 2 = 0.72) or/and sand (AOI1b, AOI3b: R 2 ≥ 0.65) contents may also significantly influence (linear correlation) the soil reflectance pattern.Spielvogel et al. (2004) [26] stated that soil spectral reflectance could be affected by (in order of priority) OM (principal determinant), CaCO3, clay, silt, and sand contents, even when no significant correlations between PCst1 and the texture parameters existed.The glaciomorphological and pedological genesis of the young morainic soil landscape of Northern Germany result in the soil reflectance not showing the darkening effect of Fe-oxides [49].Thus, for interpreting the MSPA results at the single field polygon level, additional soil texture information (e.g., soil maps) might be necessary for mapping critical interpretation zones, even at the landscape scale.
After eliminating the temporal effects, all significant linear regression models for the OM estimation were improved (R 2 ≥ 0.71) with the exception of AOI1a, AOI1d, and AOI3b.The model robustness (RMSE) varies between 8.0% and 16.5% (Table 4b).

Field Composites
The number of detected bare soil images, the number of identified soil sampling outliers, the number of remaining (used) soil samples, and the coefficient of determination (R 2 ) for each regression analysis between PCst1 and the corresponding soil parameter (OM, clay, silt, and sand) is depicted in Table 5a.Additionally, in part (b) of Table 5, the characteristics of the improved OM prediction model, such as the best percentile threshold for eliminating temporal effects, the number of sampling points used for the prediction model generation, model-quality and model-robustness criteria (R 2 and RMSE), and the results of the OM prediction model, expressed as the extremes and the range of modeled OM values are shown.
Similar to the single fields, significant (p-value <2.2 × 10 −16 ) linear (after log-transformation) regression models for PCst1 and OM with a coefficient of determination R 2 ≥ 0.6 were computed for both types of field composites.Due to the large sample sizes for the side-by-side field composites (>100 samples) and distant field composite (>500 samples) any negative effects on statistical analysis were avoided.The results reveal that no correlation between PCst1 and the corresponding soil texture parameters exists.Side-by-side field composites showed a variable R 2 (0.60 ≤ R 2 ≤ 0.77), with the lowest value for AOI3 and the highest value for AOI1.Analyzing AOI_CAL, while simultaneously testing different bare soil image combinations, demonstrated an R 2 between 0.62 and 0.67, with best result using the bare soil images from 15 September 2009, 20 September 2009, and 22 September 2010 (Table 5a).At the field polygon union step of the composite generation process, the probability of including temporal effects of vegetation and land management increases considerably due to diverse management practices (e.g., farming techniques, crop types, and schedule of arable farming).Moreover, we assume that the prediction model is also influenced by the complex cause and effect relationships of relief, soil moisture and soil texture.
For this reason, it is absolutely necessary to exclude temporal effects using the stability and percentile analysis steps of MSPA to improve the OM prediction modeling and functional soil mapping.So, all significant linear regression models for the OM estimation were improved regarding quality (side-by-side field composites: R 2 ≥ 0.63, distant field composites: 0.66 ≤ R 2 ≤ 0.69) with the exception of AOI1 Jarmen.The RMSE ranges between 7.5% and 13.8% (Table 5b).
Analogous to Blasch et al. 2015 [1], high variance between some log OM values (in the range of 0.3% and 0.7%) and PCst1 was observed, which indicates the influence of several parameters on the PCst1 spectral reflectance, even for significant regressions with high R 2 values.We assume that the complex interdependency of relief, soil texture, and soil humidity influence the soil reflectance characteristics.This might explain the variation in the prediction model quality criteria between test sites.Hypothetically, flat relief with no inclination (AOI1, AOI4), and/or homogenous soil textures (AOI2, AOI4, AOI5), have a small influence.Rather, analyzing the different physical-geographical backgrounds of all five test sites, we suspect the considerable factors are principally the heterogeneous soil texture (especially high clay content) for AOI1 Jarmen (R 2 0.77, RMSE 9.7%), the undulating/hilly relief with eroded hilltops (CaCO3 presence) for AOI2 Sassen (R 2 0.72, RMSE 9.6%), the combination of hilly relief with steep, partly eroded slopes (AOI3f: CaCO3-presence) and heterogeneous soil textures for AOI3 Buchholz (R 2 0.63, RMSE 13.8%), and the distance to the groundwater table for AOI4 (R 2 0.71, RMSE 13.7%) (Figures 2 and 3).AOI5 shows the best prediction criteria (R 2 0.81, RMSE 8.0%).For all test sites, the absolute RMSE fluctuate up to 1.3 OM-%.However, by analyzing diverse bare soil image combinations for AOI_CAL, robust improved OM prediction models with R 2 ≥ 0.66 and RMSE ≤ 7.7% were computed.The AOI_CAL test site contains all discussed physical-geographical characteristics of four MSPA-studied test sites.
Regarding prediction quality and robustness, for both single fields and field composites (except for AOI3 and AOI_CAL) the results of OM estimation based on MSPA using broadband multispectral RapidEye satellite imagery are in the same range of quality as other modeling techniques using proximal sensors.Using a narrowband VIS-NIR spectrometer for stationary and in situ proximal sensing at the field scale, Ji et al. (2015) [25] modeled soil organic carbon with R 2 = 0.71 by calibrating field spectra data against large laboratory-derived spectroscopic databases of local soils.With field spectrometry using a narrowband spectroradiometer in the VIS, NIR and SWIR range, Denis et al. (2014) [26] could improve field scale-based measurement results of soil organic carbon from R 2 0.77 to 0.89 and 0.89 to 0.91 by soil shadow corrections.In general, airborne hyperspectral imagery provide more accurate OM modeling results (e.g., R 2 ≥ 0.93 [26]).

Regional Prediction Model for OM
According to the highest R 2 = 0.69 and the very robust RMSE 7.5%, the prediction model of AOI_CAL (here 351 soil samples) with the greatest improvement based on a 60% threshold of the corresponding overall stability image, after eliminating the model input soil sampling points within the detected temporal pattern, was chosen as the best regional prediction model option for mapping the study area at the landscape scale.To produce SSM user-friendly functional soil maps based on real OM values, it is essential to re-transform the linear regression function and the coefficients of AOI_CAL back to a non-linear regression function.Consequently, the re-transformed regression equation of the regional prediction model is: y = 426245689 × x −2.094 , whereby y = OM in % and x = PCst1.Functional OM soil maps based on this regional prediction model are computed with an absolute RMSE of 1.3 OM-%.Strong tendency of underestimation exists starting from 15 OM-% on with increasing prediction error (predicted versus measured OM values).This was also confirmed at the model-independent validation test site AOI2 (absolute RMSE > 1.4 OM-%).Applying the found prediction model to all 525 soil samples of AOI_CAL and to 147 soil samples of AOI2, the prediction error varies between 1.0 OM-% (AOI_CAL) and 1.4 OM-% (AOI2) for the OM main value range of the Demmin region (0-15 OM-%).Additionally, the improved regional prediction model (Figure 5) demonstrates the accuracy of OM prediction for OM values below 15 OM-% by comparing them to measured OM values.

OM Estimation Based on the Multifield Landscape-Scale Version of MSPA
To achieve high-quality model input data for the MSPA at the landscape scale, the MSPA NDVI-threshold-adapted version of the in-house selection algorithm [47] was applied to the study area field segments shapefile.Thus, for the NDVI thresholds (MEAN 0.2, SD 0.07), three bare soil image mosaics from lowest to highest SD value were created, applying the modified algorithm to 15 suitable datasets out of the RapidEye time series (Figure 6).These bare soil image mosaics are composed of 1415 field polygons, covering an area of 343.48 km 2 .In some polygons, strong temporal patterns (e.g., vital and ripe vegetation and land management) are still visible.The NDVI-MEAN threshold controls the quantity and the quality of the bare soil image per field segment.Reducing this threshold leads to fewer field segments for bare soil image detection (e.g., for NDVI-MEAN 0.15 1151 polygons; for NDVI-MEAN 0.1 569 polygons) and to considerably-diminished strong temporal patterns, such as the presence of vital plants and land management.Despite the substantial improvement, the separation between ripe vegetation and bare soil images still remains difficult.Different spectral resolution and position of wavelength of other multispectral remote sensing sensors (e.g., Landsat-8, Sentinel-2) and hyperspectral Earth observation systems (e.g., EnMAP, HISUI, HyspIRI, PRISMA), as well as the combination of multispectral and radar imagery (e.g., ASAR/ENVISAT [50]), might be promising for improved separation of bare soil and ripe vegetation as well as for detection of land management effects.
Based on the three PCA-transformed bare soil mosaics (Figure 7) using the correlation matrix of AOI_CAL (Table 6), the enhanced regression equation of AOI_CAL was transferred to the study area fields.Thus, a functional soil map of OM was derived and classified with the adaption to the official German Soil Survey Description [42] (Figure 7).
To demonstrate the quality of the OM soil map at the landscape scale based on the regional prediction model and to visualize the behavior of the predicted OM pattern, Figure 8 shows the map sections of each test site compared with the field-scale OM soil map of each test site based on the test-site-specific prediction model, including the corresponding laboratory-analyzed OM soil sampling results.For the landscape-scale OM soil map the AOI2 test site is the primary validation site, which is not represented in the applied regional prediction model (Chapter 2.5).Consequently, the OM values of this validation test site were modeled in this map independently of the corresponding soil samples.Observing the quality criteria and robustness of the regional (R 2 0.69) and local prediction models of AOI1 (R 2 0.77), AOI2 (R 2 0.72), AOI3 (R 2 0.63), AOI4 (R 2 0.71), and AOI5 (R 2 0.81), besides the significant improvement in the model robustness (especially AOI3, AOI4), a quality decrease from the field scale to the landscape scale based OM maps was expected, with the exception of AOI3.
When checking the mapped and classified OM values, a very similar size and distribution of class zones between the soil maps were noted at the validation test site AOI2, with slight variations at AOI2 (OM classes ≥8.01%), and at the calibration test sites AOI1 and AOI5, with slight changes in the southern field areas (AOI1a, AOI1c) and strong changes in the some field areas of AOI1e and AOI1f.Generally, at AOI3, the majority of the sizes and distributions of the predicted OM class zones had only slight variations; however, in several cases, the same zones were shifted to the next higher OM class (AOI3a, AOI3b, and AOI3f).At AOI4, a strong change in the size and distribution were observed with an increase of the OM class zone from 4.01 to 8.00, particularly in the eastern field.
To assess the classification accuracy of field-scale and landscape-scale based OM maps, confusion matrixes were computed for each test site.For all test sites together, Table 7 reveals that the overall accuracy (OAA) between both map types differs about 6% (field-scale OAA: 62%; landscape-scale OAA: 56%) with similar slight tendency of overprediction in lower OM classes (≤3.00) and underprediction in higher OM classes (≥4.01).In contrast to the landscape scale classification, the field scale classification contains low underprediction for the OM class 3.01 to 4.00.The higher producer's accuracy (≥60%) indicates for both maps substantial correct class assignment of corresponding reference samples for the classes 2.01 to 8.00 at the field scale and for the classes 2.01 to 3.00, 4.01 to 8.00, and 8.01 to 15.00 at the landscape scale.Moderate producer's accuracy (40% to 60%) was obtained for the classes 4.01 to 8.00 (field-scale OM map) and 3.01 to 4.00 (landscape-scale OM map).Except for the classes ≤ 2.00 and 8.01 to 15.00, the probability that a predicted value of the classified OM map actually represents the class of measured value (consumer's accuracy) is 3% to 15% higher in the field scale OM soil map.As a measure of strength for the relationship between classification results and reference data, the calculated Kappa coefficients of agreement is 0.5 (moderate agreement) for the field scale OM map and 0.4 (moderate agreement) for the landscape scale OM map.The moderate agreement as well as the relatively low OAA result from the prediction error range (absolute RMSE: ≤1.3 OM-%) (Section 3.1.2) in relation with small class sizes of 1 OM-% for lower OM values.Considering a one-class-classification-error as acceptable, the adjusted OAA increases considerably (field-scale OAA ± 1 class: 98%; landscape-scale OAA ± 1 class: 96%).
The conventional OAAs of the validation test site AOI2 (field-scale OAA: 65%; landscape-scale OAA: 61%) as well as the calibration test sites AOI1 (field-scale OAA: 69%; landscape-scale OAA: 62%) and AOI3 (field-scale OAA: 56%; landscape-scale OAA: 50%) are 4% to 7% more accurate for the field scale maps.AOI4 shows considerable deterioration of 12% (field-scale OAA: 65%; landscape-scale OAA: 43%) due to low image quality of one input bare soil image mosaic with considerable temporal effects of ripe vegetation in the eastern part of the field at the field polygon.The OAA at AOI5 (field-scale OAA: 57%; landscape-scale OAA: 58%) improved about 1% from the field scale to the landscape scale map because of increasing correct value assignment for the class 8.01 to 15.00.

Operability of MSPA for OM Estimation Based on the Regional Prediction Model "Demmin"
According to the crop type and the image acquisition conditions (e.g., signal to noise ratio), the most suitable temporal window of RapidEye data for MSPA at the study area is generally the phase of seeding (winter crops: mid-August to October; summer crops: mid-March to mid-May).Using the found regional prediction model "Demmin" (Section 3.2), OM can be predicted at the field-scale as well as at the multi-field landscape scale for any single field or field composite of the study area Demmin without additional time-consuming, cost-intensive soil sampling and laboratory analyses.
First, exactly three bare soil images of good quality without any vegetation and/or land management influence are needed.Depending on site-specific characteristics and land management practices, images should not be more than three to five years apart in terms of time, to avoid any OM degradation effects.To create a precise multitemporal bare soil image of the agrarian field of interest without any mixed-pixel influence of field boundary vegetation and/or surrounding trees, all three bare soil images should be stacked together and clipped to the precise extent of the field boundary with a buffer of 20 m.The correlation matrix of AOI_CAL (Table 6) must be applied as transformation parameters to the standardized PCA of the multitemporal stacked bare soil image, to obtain the field-specific PCst1-image as prediction model input data.Finally, OM soil map can be produced by introducing the PCst1-image to the regression equation of the regional prediction model (Section 3.2) using for example the Raster Calculator of ArcGIS Spatial Analyst.
Due to the commercial character of the RapidEye satellite system (0.95 €/km 2 ) [51], the 54 used RapidEye datasets would cost 42,784.20€ for the study area of 834 km 2 .Covering an area of 343.48 km 2 , the generated OM soil map at the landscape scale based on the regional prediction model "Demmin" cost 124.56€/km 2 (namely, 1.25 €/ha) exclusive of additional costs of manpower, working time, soil sampling, and laboratory analyses.In case of the regional model application to any field of the study area, soil sampling, and corresponding laboratory analyses would be redundant.Assuming 20 RapidEye datasets (ordered within the crop type-depending, most suitable temporal window) as sufficient for MSPA application, the costs would range between 4.75 €/ha and 0.95 €/ha for common farm sizes (2000 ha to 10,000 ha) in Mecklenburg-Western Pomerania, taking into account the minimum order size of RapidEye data (500 km 2 ).
In the future, the combination of the MSPA method with free available Sentinel-2 or Landsat-8 data will enable a high operational use at multiple scales and no cost for image acquisition.This could be improved with the link to locally-developed and/or existing soil spectral libraries such as the European Land Use/Cover Area Statistical Survey database (LUCAS) [52].

Conclusions and Outlook
This research presents the development and implementation of a multi-field landscape-scale version of Blasch's et al. (2015) [1] field-scale based multitemporal soil (reflectance) pattern analysis (MSPA) method for organic matter (OM) prediction using RapidEye satellite imagery.In the process of method development, the field-scale based MSPA was successfully applied to single fields and field composites with diverse physical-geographical location characteristics.Subsequently, functional soil maps based on a field-specific non-linear relationship between OM soil sampling data and the soil pattern detecting first principal component (PCst1) could be computed.By combining the field-specific location characteristics of single fields and field composites, a regionally significant and very robust OM prediction model (R 2 0.69; relative RMSE 7.5%) with a prediction accuracy of 1.3 OM-% (absolute RMSE) based on 351 soil samples within highly-stable spatiotemporal reflectance patterns was found for the study area Demmin.The transfer of the regional OM prediction model to the multi-field landscape scale was realized using the extracted correlation matrix from the standardized principal component analysis step and the bare soil image mosaics generated by the modified selection algorithm [47].This enabled the generation of a precise OM soil map at the landscape scale for all agrarian fields within the study area.Resulting landscape-scale based OM soil map demonstrated a high degree of accordance to field-specific OM soil maps.
Consequentially, the MSPA approach using the found regression equation of the regional prediction model "Demmin" combined with RapidEye data allow precise OM estimation for agrarian fields of the Demmin region, as cost-efficient, highly operable and rapid method.All application steps as well as the most suitable temporal window of RapidEye data have been explained in detail.To evaluate the validity of the regional prediction model "Demmin" for distant agrarian fields of the young morainic soil landscape of Northeastern Germany, further studies with additional soil sampling will be conducted on agrarian fields of the Quillow catchment area, which is also located in the northeastern lowlands of Germany, 5 km west of Prenzlau and ~75 km southeast of Demmin.
As Blasch et al. (2015) [1] already stated, in the future, high-resolution functional OM soil maps at multiple scales that are derived from static soil patterns could be applied in the form of ancillary data for SSM in PA (e.g., sowing intensity) and for further soil science related activities (e.g., soil survey preparation, digital soil mapping, erosion models, and organic carbon stock modeling).

Figure 1 .
Figure 1.The digital elevation model (resolution 10 m) [29] of the study area (red rectangle), study area field segments (grey polygons) and the five test sites (AOI 1-5) composed of 16 agrarian fields (black polygons).

Figure 2 .
Figure 2. Physical-geographical characteristics of each single field test site with soil sample locations: (a) RapidEye bare soil image; (b) DEM; (c) Soil texture map).

Figure 3 .
Figure 3. Physical-geographical characteristics of each field composite test site with soil sample locations: (a) RapidEye bare soil image; (b) DEM; (c) Soil texture map.

Figure 4 .
Figure 4. Analyzing steps of the MSPA workflow.

Figure 5 .
Figure 5. (a) Regional prediction model-Scatterplot of the relationships between OM and the standardized PC 1 of the multitemporal bare soil images of AOI_CAL based on the static soil pattern (threshold: 60%); (b) predicted versus measured OM values for OM values below 15 OM-% based on the regression equation of the regional prediction model (brown line: One-to-one line; orange line: absolute RMSE boundary).

Figure 6 .
Figure 6.Bare soil mosaics for all fields of the study area (red rectangle) including AOI_CAL (yellow polygons/circle) and the validation test site AOI2 (red polygons/circle).

Figure 7 .
Figure 7. (Left) PCst1-transformed bare soil image mosaics of the study area (red rectangle) using the correlation matrix of AOI_CAL (yellow polygons/circle); (right) OM functional soil map based on the regional regression equation (Chapter 3.2) and the PCst1-transformed bare soil image mosaics (OM classification adapted to the official German Soil Survey Description with soil sampling results); validation test site AOI2 (red polygons/circle).

Figure 8 .
Figure 8. (Left) Map section of the landscape-scale OM soil map based on the regional prediction model for each test site with the validation test site AOI2 (red rectangle); (right) field-scale OM soil map of each test site based on the test-site-specific prediction model (OM classification adapted to the official German Soil Survey Description with soil sampling results).

Table 1 .
Baseline data and specific data characteristics.

Table 2 .
Soil sampling design for each single field and field composite test site.

Table 3 .
Descriptive statistics of soil parameters (SD: standard deviation; CV: coefficient of variation).

Table 4 .
(a) Results of the identified bare soil images, outlier detection and R 2 for the relationship between the soil parameter (dependent variable) and the PCst1 of the multitemporal images (independent variable) for each single field test site and each single field polygon of the field composite test sites; (b) characteristics and results of improved OM model based on MSPA; (* sample number for texture analysis; B.S.I = Bare soil images; P. = percentile).

Table 5 .
(a) Results of the identified bare soil images, outlier detection and R 2 for the relationship between the soil parameter (dependent variable) and the PCst1 of the multitemporal images (independent variable) for each field composite (site-by-site: AOI1, AOI3; distant-field: AOI_CAL); (b) characteristics and results of improved OM model based on MSPA; (* sample number for texture analysis; B.S.I = Bare soil images; P. = percentile).

Table 7 .
Confusion matrix of the measured OM values (reference) with predicted OM values (classification) for field-scale and landscape-scale based OM soil maps (P.A.: producer's accuracy, C.A.: consumer's accuracy).