Next Article in Journal
Mango Fruit Processing: Options for Small-Scale Processors in Developing Countries
Next Article in Special Issue
Estimation of Soil Nutrient Content Using Hyperspectral Data
Previous Article in Journal
The Joint Action of Some Broadleaf Herbicides on Potato (Solanum tuberosum L.) Weeds and Photosynthetic Performance of Potato
Previous Article in Special Issue
Evaluation of Deep Learning for Automatic Multi-View Face Detection in Cattle
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Agriculture
Volume 11
Issue 11
10.3390/agriculture11111104
agriculture-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Simplified and Hybrid Remote Sensing-Based Delineation of Management Zones for Nitrogen Variable Rate Application in Wheat

by
Mohammad Rokhafrouz
1,
Hooman Latifi
1,2,*,
Ali A. Abkar
3,
Tomasz Wojciechowski
4,
Mirosław Czechlowski
4,
Ali Sadeghi Naieni
5,
Yasser Maghsoudi
1 and
Gniewko Niedbała
4
1
Department of Photogrammetry and Remote Sensing, Faculty of Geodesy and Geomatics Engineering, K. N. Toosi University of Technology, Tehran 15433-19967, Iran
2
Department of Remote Sensing, University of Würzburg, Oswald Külpe Weg 86, 97074 Würzburg, Germany
3
AgriWatch BV, World Trade Center Twente, Industrieplein 2, 7553 LL Hengelo, The Netherlands
4
Department of Biosystems Engineering, Faculty of Environmental and Mechanical Engineering, Poznań University of Life Sciences, Wojska Polskiego 50, 60-627 Poznań, Poland
5
Photogrammetry and Remote Sensing Department, National Cartographic Center, Tehran 1684-13185, Iran
*
Author to whom correspondence should be addressed.
Agriculture 2021, 11(11), 1104; https://0-doi-org.brum.beds.ac.uk/10.3390/agriculture11111104
Submission received: 20 September 2021 / Revised: 26 October 2021 / Accepted: 2 November 2021 / Published: 5 November 2021
(This article belongs to the Special Issue Digital Innovations in Agriculture)
Browse Figures
Versions Notes

Abstract

:
Enhancing digital and precision agriculture is currently inevitable to overcome the economic and environmental challenges of the agriculture in the 21st century. The purpose of this study was to generate and compare management zones (MZ) based on the Sentinel-2 satellite data for variable rate application of mineral nitrogen in wheat production, calculated using different remote sensing (RS)-based models under varied soil, yield and crop data availability. Three models were applied, including (1) a modified “RS- and threshold-based clustering”, (2) a “hybrid-based, unsupervised clustering”, in which data from different sources were combined for MZ delineation, and (3) a “RS-based, unsupervised clustering”. Various data processing methods including machine learning were used in the model development. Statistical tests such as the Paired Sample T-test, Kruskal–Wallis H-test and Wilcoxon signed-rank test were applied to evaluate the final delineated MZ maps. Additionally, a procedure for improving models based on information about phenological phases and the occurrence of agricultural drought was implemented. The results showed that information on agronomy and climate enables improving and optimizing MZ delineation. The integration of prior knowledge on new climate conditions (drought) in image selection was tested for effective use of the models. Lack of this information led to the infeasibility of obtaining optimal results. Models that solely rely on remote sensing information are comparatively less expensive than hybrid models. Additionally, remote sensing-based models enable delineating MZ for fertilizer recommendations that are temporally closer to fertilization times.

1. Introduction

In recent years, there has been an intense growth in the world population, which is projected to reach 9.7 billion by 2050 [1]. Population growth puts enormous pressure on agricultural productivity growth, but also on the increasing environmental impact of the agri-food sector [2,3]. In many regions of the world, small farms are the main food producers, and this group will be under pressure to increase production efficiency [4]. Important elements in the process of agricultural production growth are biological and technological progress [5]. In the scope of technological progress, the most promising developments are considered to be precision agriculture (PA) (agriculture 3.0) and digital agriculture (agriculture 4.0) [6,7,8].
Precision agriculture or site-specific management can provide food security and sustainable development [3,6,9,10,11]. It is based on innovative system approaches which comprise several technologies, such as global navigation satellite system (GNSS), geographic information system (GIS), proximal sensing (PS) and remote sensing (RS), artificial intelligence (AI), machine learning (ML), automatic guidance, section control, variable rate technology (VRT) and advanced information processing for timely within- and between-season crop managements [12,13,14,15]. The main purpose of PA is to optimize crop management concerning spatial and temporal variabilities, which results in optimized utilization of farm inputs such as fertilizers, pesticides, herbicides and seeds [16]. All of this is aimed at increasing farm profitability and achieving the Sustainable Development Goals (SDGs) such as No Poverty, Zero Hunger and Reduced Inequalities. For such a purpose, a wide range of data and information from field inventory, crop growth and yield patterns must be analysed [17]. With this correctly processed information, agricultural inputs such as fertilizers, water or energy can be applied in a spatially variable manner using homogeneous production zones, i.e., management zones (MZ) [16].
The RS and PS are the most researched technology in PA [18]. Their effective use requires the delivery of information and communication technologies (ICT) tools, including algorithms for RS, that are useful to the user: farmers with limited knowledge capital, which is currently indicated as one of the main barriers to the adaptation and spread of PA [18]. The need to support and develop simplified, low-tech precision farming methods seems therefore justified.
MZ are defined as sub-units of farm fields with a relatively homogeneous combination of yield-limiting factors [16,19]. Each zone can be managed with a different but specific single-rate management practice to maximize the efficiency of farm inputs [16,19]. Methods to create MZ have been developed for almost 3 decades, and the evolution of methods to create them is widely reported in the literature [2]. Generally, MZ delineation approaches can be categorized based on the provided data and information from different sources [2]. These methods are generally based on farmers’ knowledge [20], soil physical and chemical attributes [21,22,23,24], geomorphology [25], yield [26,27,28,29,30,31], electrical conductivity (EC) [32,33] and RS [17,34,35] data, and also hybrid models that combine information from different data sources [36,37,38,39,40,41,42,43,44,45,46,47].
Each of those approaches has its pros and cons. Although the hybrid method is theoretically comprehensive and more accurate, additional field measurements such as soil sampling or proximal soil sensing are expensive, labor-intensive and time-consuming, and require seasonal sampling to specify nutrient level due to temporal variability of soil properties [22,48,49,50]. Besides, large commercial agricultural fields reportedly do not completely represent spatial variability [51]. RS methods deliver key components of precision farming and provide valuable data from crop coverage and actual crop growth patterns to delineate MZ [52]. With spatiotemporal continuity along with cost-effectiveness, RS has a capacity for time-series analysis [17]. However, optical satellite imagery is associated with the main limitation of being affected by atmospheric haze or cloud, which often occurs in temperate and rainy regions [17].
A broad range of active/passive satellite RS data is currently available with various properties, such as spatial resolution, temporal resolution, spectral range and viewing geometry [53]. The launch of Sentinel-2A (2015) and Sentinel-2B (2017) satellites by the European Space Agency (ESA) boosted the PA applications since the data are freely available [54]. Sentinel-2 satellites are equipped with a passive multispectral instrument (MSI), including 13 spectral channels, 4 bands at 10 m, 6 bands at 20 m and 3 bands at 60 m spatial resolution [55]. These satellites have a high revisit time of ten days with one satellite and five days as a constellation (2–3 days at mid-latitudes) [56]. Currently, ESA provides the Level-2A products of this mission as bottom of atmosphere (BOA) reflectance and atmospherically corrected images [57].
Remote sensing of vegetation is mainly based on the green (495–570 nm), red (620–750 nm), red-edge (680–730 nm), near- and mid-infrared bands (850–1700 nm) regions of the spectrum [35]. In order to obtain information from vegetation status, these bands can be used to derive vegetation indices [35]. The normalized difference vegetation index (NDVI) [58], a normalized difference between the reflectance of red and near-infrared (NIR) spectral bands, is the most common crop parameter used in MZ delineation, because of its ease of calculation and interpretation [2,35]. Along with NDVI, some biophysical variables such as leaf area index (LAI), fraction of absorbed photosynthetically active radiation (FAPAR) and the fraction of vegetation cover (FVC) were also used in this research [59]. LAI is defined as half the developed area of photosynthetically active elements of the vegetation per unit of horizontal ground area [59], FAPAR corresponds to the fraction of photosynthetically active radiation absorbed by the canopy [59] and FVC is the ratio of the vertically projected area of vegetation to the total surface area [60].
Several studies for MZ delineation on multispectral satellite images exist in the literature. Song et al. [37] delineated and compared MZ based on soil data, yield data, and crop RS information from one multispectral satellite scene, as well as their combination. Additionally, Martins et al. [47] generated a MZ map with a combination of soil attributes, EC, yield maps and a vegetation index (VI) of one multispectral satellite image. Georgi et al. [17] developed a segmentation algorithm for generating MZ from within-field crop patterns using only multi-temporal, multi-spectral satellite images.
The above examples [17,37,47] reveal that several studies conducted MZ delineation with various methods and different data types. However, none of them considered agronomy (e.g., BBCH stage) and especially climate information (e.g., soil moisture conditions, drought), which could negatively affect the implemented models. Additionally, comparative approaches need to be tested to improve the MZ delineation considering data and/or knowledge availability, time and cost-benefit analysis and accuracy.
Above all, the extent of available agricultural spatial data for farmers is also an important issue in the development of homogeneous MZ and homogeneous productivity zones. The ideal scenario is when the producer/farmer has multi-years spatial data on yield, crop vegetation, soil and climate, but literature reports indicate that the level of farm datafication varies [61]. Practice shows that the ideal scenario exists in regions with high adaptation of PA technologies, e.g., where the use of combine harvester yield monitors is common practice [18]. In most cases, these data are less available. This also applies to smaller farms with lower adaptation of PA technologies. Proposing new digital technologies for such users, such as Decision Support Systems (DSS) and Farm Management Information Systems (FMIS), there is a need to develop simplified algorithms/methods for the creation of MZ, based on one or two seasons’ data. The problem of determining the MZ of simplified algorithms that work in practice, i.e., in FMIS, is presented in research by Santaga et al. [62].
The overall aim of the study is to generate and compare MZ maps prepared using different models. The first model was termed as “RS- and threshold-based clustering” and was adopted from Georgi et al. [17] with partial modifications. The second model was a “hybrid-based and unsupervised clustering” model, a hybrid model in which a combination of data from different sources was utilized for MZ delineation. The last approach was called “RS-based unsupervised clustering”, which was similar to the first model but with a different classification. The secondary objectives were to statistically analyse the models’ outputs, as well as to improve the accuracy of MZ delineation by incorporating agronomy and climate information in the applied models. Finally, a MZ map was presented to guide spring mineral nitrogen fertilization of winter wheat in specific phenological phases based on fertilization dates.

2. Materials and Methods

2.1. Study Area

The research was conducted on the 50.2 ha field on the experimental farm of the Poznań University of Life Sciences, RGD Brody, located in Brody (52.43 N, 16.29 E according to WGS84), Wielkopolskie Voivodeship, Poland (Figure 1). The field dedicated to the research was covered with winter wheat crop (Triticum aestivum cv. ‘RGT Reform’), carried out in the reduced soil tillage system and non-irrigated.
The average annual precipitation sum of the study area (1960–2019) was 599 mm, and the annual mean air temperature was 8.5 °C, while in 2019 and 2020, the annual precipitation sum was 462 and 520 mm respectively, with the mean air temperature of the study at 10.8 and 10.6 °C, respectively. Meteorological data were obtained from the meteorological station of Brody Experimental Station and were recorded according to the World Meteorological Organization guidelines.
The farm soils in the study area are light, loamy sands, developed on loamy sands overlying loamy material, and are classified as Albic Luvisols according to World Reference Base nomenclature [63,64].

2.2. Data

2.2.1. Soil Sampling

Soil physical and chemical properties were determined by infield sampling and laboratory analysis. Soil samples were taken mechanically, in a semi-automatic operation, from the 0–30 cm field layer in a 4 ha grid, where one average, mixed sample consisted of 16 primary samples. Sampling was carried out in spring 2020. In the present study, analyses were performed to determine pHKCl, phosphorous (P2O5), potassium (K2O) and magnesium (Mg) contents. These are the most commonly used parameters determined by farmers in agrotechnical practice, especially small-scale farming, because of the low implementation costs. Additionally, an analysis of soil organic matter (OM) content was performed as one of the important factors determining sorption and water properties of arable soils.
Soil physico-chemical parameters were determined as follows: pH in 1M KCl according to PN-ISO 10390-1997, available phosphorus and potassium by the Egner-Riehm method [65] according to PN-R-04023 and PN-R-04022 respectively, and magnesium by the Schachtschabel method [66] according to PN-R-04020. OM content was determined by the Tiurin method (T methode) [67] with the Van Bemmelen coefficient of 1.724.

2.2.2. Yield Data

Winter wheat yield data were recorded automatically during harvest in 2019 and 2020 with a modified Claas Lexion 480 combine harvester. Data were recorded at a temporal resolution of 1 Hz for each of the harvester passes. The recorded raw yield data were post-processed and filtered to mitigate lag times and exclude outliers. A detailed description of the combine harvester prototype equipped with a system for monitoring qualitative and quantitative grain parameters is described by Czechlowski and Wojciechowski [68,69].

2.2.3. Elevation Data

The basic hypsometric data were obtained using the measuring system of the modified combine harvester described by Czechlowski and Wojciechowski [64]. The combine was equipped with a Novatel RT2 PROPAK V3 GNSS receiver with a GPS-702-GG: a dual-frequency (L1/L2) antenna and a SmallTRIP 3.2 GPRS/NTRIP modem with automatic connection to the Real-Time Kinematic (RTK) NAWGEO service of the ASG-Eupos network. The possibility of using this type of data to create digital elevation models of agricultural field surfaces was reported by Czechlowski et al. [70]. The slope map was generated based on the interpolated elevation map (DEM) as a percent slope (See Preprocessing Section).

2.2.4. RS Data

A time series of Sentinel-2 L2A images from 1 January 2018 to 1 July 2020 was downloaded from The Copernicus Open Access Hub [71]. Sentinel-2 L2A data are atmospherically and geometrically corrected. Additionally, layers such as a scene classification layer, cloud mask, cloud shadow mask and snow mask were provided along with Sentinel-2 L2A raw data. During the downloading stage, it was attempted to download the data with minimum possible cloud cover by checking the quick layer of every data point. The quick layer of each data point is accessible in the Copernicus Open Access Hub which is embedded in the detailed information of each Sentinel-2 scene. Finally, 119 scenes of Sentinel-2 L2A were downloaded. Table 1 shows the number of downloaded data per year and month.

2.3. Models

Three models were separately implemented to delineate MZ. In the following sections, each model is described in detail. First, a brief introduction of each model is provided, then its flowchart is presented and related parts are explained.

2.3.1. Model-1 (RS- and Threshold-Based Clustering)

The first model was fundamentally adopted from Georgi et al. [17], who developed a segmentation algorithm for generating MZ of within-field crop patterns by solely using multi-temporal and multi-spectral satellite images. Thus, the input to the model was the time series of RS data. However, the RapidEye data used in [17] was replaced with Sentinel-2 L2A data because of its free and open data policy. Moreover, a different approach was applied for selecting cloud-free and cloud shadow-free data using mask layers embedded in Sentinel-2 L2A products. The workflow of model-1 is summarized in Figure 2 and the whole process of this model was subdivided into 4 parts (Sentinel-2 data processing, Data selection, Processing of NIR bands, and Segmentation and classification), and a detailed description of each part is provided in the following sections.

Sentinel-2 Data Processing

As inputs to the model, 119 scenes of the Sentinel-2 L2A were imported. First, all bands were resampled to 10 m pixel size, which was inevitable due to differences in spatial resolutions. The cloud-free and cloud shadow-free data were extracted by investigating the quality scene classification layer, quality cloud confidence, cloud probability mask, cloud mask and shadow mask layers embedded in Sentinel-2 L2A products. For each scene, the mentioned layers were investigated by visual quality control. Finally, NDVI was calculated for all images. The entire workflow was conducted in SNAP V7 software [72].

Data Selection

The data were selected with two constraints. Standard deviation (SD) was calculated for each NDVI data point, and the SD values <0.02 were dropped to exclude the images of dense vegetation cover with no spatial patterns. Then, the mean of each NDVI data point was calculated, and the values between 0.3 and 0.78 were selected since the values <0.3 and >0.78 depict the corresponding images of vegetation canopy, which is too sparse (bare soil background) or too dense, respectively [73]. When the canopy becomes too dense, NDVI saturates because red reflectance does not change much, but near-IR reflectance increases [73]. This stage was implemented in Python 3.8 by using ‘rasterio’, ‘unidip’ and ‘numpy’ packages. For selected dates, NIR bands of Sentinel-2 data were extracted, on which the consequent steps were implemented. This model was conducted based on NIR bands, while NDVI data were used only for data selection since ratio indices such as NDVI cause noise patterns and artifacts that challenge the MZ delineation [17,74]. This resulted in 24 out of 119 raster data, as summarized in Table 2.

Processing of NIR Bands

Each image was converted to relative values by Equation (1), i.e., a normalization to a percentage, where 100% was equal to the average NIR value of each image.
N o r m e d   p i x e l   v a l u e = ( P i x e l   v a l u e − M i n i m u m   M e a n − M i n i m u m   ) × 100
where Minimum is the minimum value of the whole scene and Mean is the mean value of the whole scene. Then, an average of NIR bands for each year was calculated, thus the NIR time series of each year formed a raster stack. Then, normalization and averaging for the generated data were applied as well.

Segmentation and Classification

A 3 × 3 median filter was applied to eliminate the small zones and smooth the class boundaries, and the result was normalized as previously mentioned. Eventually, a thresholding method was implemented on the classification, whereby the 10%, 35%, 65% and 90% quantiles were calculated, and the final raster was classified into five classes. These quantile values were empirically chosen [17]. The final five classes were termed ‘very low’ (1), ‘low’ (2), ‘average’ (3), ‘high’ (4) and ‘very high’ (5), which correspond to yield expectancy. The processing of NIR bands, classification and mapping was conducted in QGIS V3.10 [75].

2.3.2. Model-2 (Hybrid-Based, Unsupervised Clustering)

The Second approach is based on a hybrid model for MZ delineation that combined data from different sources (see Figure 3). This method was recently applied by researchers in several studies [44,45,46,47]. However, in these studies [44,45,46,47], principal component analysis (PCA) was utilized to reduce data dimensionality and minimize the dependencies among variables. In this study, PCA was replaced with machine learning-based feature selection (random forest (RF) feature importance). The flowchart of model-2 is shown in Figure 3 and the whole workflow of this model was subdivided into 5 parts (Input data, Preprocessing, Processing, Output (MZ map) and Validation), and a detailed description of each part is provided in the following sections.

Input Data

For this method, data from 4 different sources were integrated, including (1) soil nutrition data comprising soil pHKCl, P2O5, K2O, Mg and OM, (2) topographical data comprising elevation and slope of the study area, (3) yield data from 2019 and 2020 and (4) RS data comprising NDVI and biophysical variables of LAI, FAPAR and FVC of Sentinel-2 data at the heading stage of wheat (15 May 2020 for this study area) [37]. The NDVI and biophysical variables were derived in SNAP V7 [72].

Preprocessing

Descriptive statistics (minimum, maximum, mean, SD, standard error (SE), coefficient of variation (CV), skewness and kurtosis) of soil, elevation and yield samples were calculated. Since the locations of soil, yield and elevation data are different, no geostatistical analysis was possible prior to correlation analysis and feature selection. Thus, semi-variogram parameters (nugget (C0), sill (C + C0) and range) were estimated to represent the spatial distribution of soil, yield and elevation data [76]. Several semi-variogram models, including circular, spherical, tetraspherical, pentaspherical, exponential, Gaussian, rational quadratic, hole effect, k-Bessel, J-Bessel and stable, were evaluated. The best-fit model with the lowest root-mean-square (RMS) error was selected for each data point. Then, the data were interpolated using the best-fit model and ordinary kriging (OK) procedure. Since the resolution of the satellite data is 10 m, all the variables (e.g., soil, yield and elevation data) were interpolated to 10 m spatial resolution. Based on Martins et al. [44], variables that are temporally stable and correlated with crop yield were selected, which are crucially significant to delineate MZ [47]. Therefore, a correlation matrix was generated using Spearman’s correlation to specify the relationship among variables on interpolated data with the spatial resolution of 10 m. Unlike other recent studies [44,45,46,47], in this study, RF with variance reduction criterion was utilized instead of PCA to rank the features, reduce data dimensionality and minimize the dependencies amongst variables. The final variables were selected by considering Spearman’s correlation matrix (using correlation coefficient criterion) and RF feature importance (using mean squared error criterion). Deriving descriptive statistics, correlation analysis and feature selection was performed in Python 3.8 using ‘pandas’, ‘scipy’, ‘rasterio’, ‘numpy’, ‘matplotlib’, ‘seaborn’ and ‘sklearn’ packages. The geostatistical analysis and interpolation were conducted in ArcGIS V10.7 [77].

Processing

MZ delineation was performed by the fuzzy c-means algorithm using Management Zone Analyst (MZA) software V1.0. [36,78]. Furthermore, two types of cluster validity functions, fuzzy performance index (FPI) [79,80] and normalized classification entropy (NCE) [81], were used to determine the optimum number of zones. The FPI is a measure of the degree of separation, i.e., fuzziness, between classes, with values ranging from 0 to 1 [36]. Additionally, the NCE measures the degree of disorganization between classes [36]. The minimum values of these indices suggest the optimum number of clusters since it represents the least membership sharing (FPI) or the highest amount of organization (NCE), as shown in Equations (2) and (3) [36]. The settings used in the MZA software included similarity measure = Mahalanobis distance, fuzziness exponent = 1.3, the maximum number of iterations = 300, convergence criteria = 0.0001, minimum number of zones = 2 and maximum number of zones = 8. Finally, the ideal number of zones was selected by considering the lowest values of FPI and NCE:
F P I = 1 − c ( c − 1 ) [ 1 − 1 n ∑ k = 1 n ∑ i = 1 c ( u i k ) 2 ]
N C E = c ( n − c ) [ − 1 n ∑ k = 1 n ∑ i = 1 c u i k l o g a ( u i k ) ]
where c is the number of clusters, n is the number of observations, uik is the fuzzy membership and loga is the natural logarithm. Following clustering, mapping was conducted in QGIS V3.10 [75].

2.3.3. Model-3 (RS-Based, Unsupervised Clustering)

The approach applied here differed from model-1, in that a K-means clustering algorithm was conducted to compare the classification of this model (threshold-based clustering) with a simple clustering procedure. Other components of model-3 were similar to model-1 (Figure 2). As can be seen from Figure 4, the whole workflow of model-3 was subdivided into 3 sections (Sentinel-2 data processing, Data selection and Classification) and the description of each part was provided in Section 2.3.1.
The K-means algorithm was performed with 5 classes and 100 iterations in SNAP V7 [72]. Moreover, mapping was conducted in QGIS V3.10 [75].

2.4. Model Improvement

To improve the result of the final MZ maps, the RS-based models (model-1 and model-3) were enriched with climate and agronomy information. The agronomic information used is the date of sowing and the dates of subsequent spring mineral nitrogen fertilization treatments, which are carried out three times, as is typical for this region. The timing of the three application rates is linked to the phenological phases of plant development and to climate and soil conditions (beginning of vegetation, stem-shooting phase, earing). Agronomic and drought information were used to select RS images. This is a knowledge-based selection of input RS data under new climate conditions, such as drought. Additionally, both RS-based models were performed concerning phenological phases, which starts from seeding to harvesting time using expert knowledge. To overcome this, the models should be run with single-year RS data and considering phenological phases in a specific year. The abnormal climate condition, such as drought in the analysed period, can impact the final yield map, which affects the performance of the models. The RS-based models were performed only for 2020 RS data to avoid this, by selecting data after 19 September 2019 (seeding date) in the 2019–2020 season. This resulted in 5 RS datasets that passed the data selection constraints (Table 2). Since this analysis was conducted for single-year data, the generated MZ maps can guide fertilization before the fertilization time approaches considering the growth stages of wheat. In this study area, fertilization was performed at four dates in the 2019–2020 season (23 January 2020, 17 February 2020, 25 March 2020 and 21 April 2020). With regard to fertilization dates and available selected RS data, two MZ maps for fertilization were generated for 25 March 2020 and 21 April 2020.

2.5. Sampling for Validation

A stratified random sampling procedure was conducted for drawing validation data. Yield maps were converted to relative values (see the Processing of NIR Bands Section) and then averaged based on the models over the available years [17]. For each sample, the relative yield value and corresponding class ID were sampled based on the MZ map [17]. The size and conjectured SD of each class were considered to determine the sample size of each class. The number of samples was computed by Equation (4) [82]:
N = ( ∑ i = 1 W i S i S o ) 2
where N is the number of samples, Wi is the proportion of mapped area for class i, Si is the SD of stratum i, So is the expected SE of overall accuracy and C is the total number of classes. The Si of each class was conjectured since no specific data for samples were available to determine the Si of each class by the assumption that Si is higher for classes with low area proportion. Afterwards, So was assumed equal to 0.01 [82]. Finally, equal distribution (EDi) (Equation (5)) and weighted distribution (WDi) (Equation (6)) of each class’ samples were calculated to determine the sample size of each class. The final number of samples for each class (Ni) (Equation (7)) was assessed by averaging these distributions. The sampling procedure was performed in QGIS V3.10 [75].
E D i = N × W i
W D i = N C
N i = ( E D i + W D i ) 2

2.6. Validation

To validate and evaluate the final MZ, a variety of statistical tests were performed. The Paired Sample T-test, Kruskal–Wallis H-test and Wilcoxon signed-rank test were applied to samples of each class and compared with each other to explore whether there were statistically significant differences. The Paired Sample T-test requires normally distributed data, so it was performed on logarithmically transformed sample values for the second time after the test was performed without normalization. The purpose of MZ delineation is to classify the wheat parcel into homogeneous zones, thus the separability of final zones was tested. Finally, boxplots for each model were used along with fitted lines through the medians of each class. The validation was conducted in Python 3.8 by using ‘pandas’, ‘scipy’, ‘numpy’ and ‘matplotlib’ packages.

3. Results

3.1. Model-1

Figure 5a shows the delineated, 5-class MZ map based on model-1. The higher the number of zones, the better the crop vitality and yield expectancy will be. Based on this map, the north and south-west parts of the field (Zones 5 and 4) show a more productive crop pattern compared with the west and east parts. The statistical tests (Table S1 see the (Supplementary Material)) indicated inseparability between classes for values with p > 0.05, and these values are highlighted in red in Table S1. According to Table S1, although the result of the Kruskal–Wallis H-Test, which compares the separabilities of the zones all at once, showed that all classes are separable with the p-value of 2.9 × 10−54 in all other tests, the pairs of the zones 3–4, 3–5 and 4–5 are inseparable. The results of the Paired Samples T-test for the pairs of the zones 3–4, 3–5 and 4–5 were 0.60, 0.63 and 0.95, respectively. The p-values of the Paired Samples T-test (log of data) for the mentioned zones were equal to 0.86, 0.80 and 0.70, respectively. Finally, the results of the Wilcoxon signed-rank test for these zones were 0.34, 0.08 and 0.77. Moreover, in the Wilcoxon signed-rank test, the pair of the zones 1–2 did not support the separability hypothesis, with a p-value of 0.08. On the other hand, values with p < 0.05 indicate separable zones. Additionally, the boxplot shown in Figure 6a confirms the result of the statistical tests, with overlaps observed for the pairs of the zones 3–5 and 1–2.

3.2. Model-2

The results of descriptive statistics for soil, yield and elevation samples are summarized in Table 3. Despite the sufficient number of samples for elevation (DEM), OM and yields for 2020 and 2019, the numbers of soil pHKCl, P2O5, K2O and Mg samples were fewer than the expected number for interpolation as there were just 14 samples in the 50 ha area. The authors of this study are well aware of the small sample size, but this is the typical soil sampling density used in state public advisory practice. As shown in Table 3, the average yield in 2020 dropped by approximately 400 kg ha−1 compared with that of 2019. This lower productivity is attributed to the severe spring drought that occurred in major parts of Poland, including this study area [83]. The CV of the analysed attributes can be categorized from low (CV < 12%) to moderate (12% ≤ CV < 60%) based on the classification suggested by Warrick and Nielsen [84].
Table 4 represents the results of the geostatistical analysis. It suggested the best-fit models to be exponential (DEM, yield 2019, yield 2020) J-Bessel (OM, pHKCl), Gaussian (K2O, Mg) and Hole Effect (P2O5), based on the minimum RMSE. The values of Nugget/Sill could be used to determine the degree of the spatial autocorrelation, in which the values < 25%, 25–75% and >75% suggest strong, moderate and weak spatial dependencies, respectively [85]. DEM and OM showed strong spatial dependence, while other parameters showed moderate degrees. The range of the semi-variogram was the distance over which the samples are correlated with each other [86]. A low value of Nugget/Sill and a high range of an attribute generally indicate that high precision can be obtained by kriging [85].
Figure 7 shows the maps using the best-fit model, including interpolated soil, elevation and yield, along with RS data.
The maps suggested a high correlation of RS data with the yield map of 2020. Moreover, high values were generally observed in the central part of the field. Soil pHKCl and P2O5 maps were consistent since the values of P2O5 were high and neutral in the west and east parts of the field. However, the values were low and somewhat acidic in the southern part. In terms of K2O, Mg and OM, higher values were observed in the Eastern part. Besides, DEM and slope values were higher in the western part. Figure 8 and Figure 9 show the results of correlation analysis and selection of features with high correlation with yield 2020 data. The RS data (NDVI, LAI, FCV and FAPAR) were highly and positively correlated with yield 2020. The feature selection also showed that yield 2019 and NDVI are appropriate features for clustering.
Besides, the optimum number of classes was found by computing two cluster validity indices (FPI, NCE). Figure 10 shows the plotted values of FPI and NCE against the number of clusters, with the optimum number being the value at which FPI and NCE are minimum, i.e., 5 clusters. This was consistent with model-1 results in terms of the number of zones.
The result of MZ in 5 classes by fuzzy c-means clustering (Figure 5b) showed a better condition in the central part of the field from west to east (Zone 5). This map was consistent with the features (NDVI and yield 2019 data) that were used in the MZ delineation process. As reported in Table S2, the p-value of the Kruskal–Wallis H-Test for this model was 3.2 × 10−117, which means that all of the zones support the separability hypothesis. However, the results of other statistical tests showed that the pair of the zones 2–3 did not support the separability hypothesis, with p-values of 0.37 (Paired Samples T-test), 0.46 (Paired Samples T-test (log of data)) and 0.55 (Wilcoxon signed-rank test). This can also be observed in the boxplot (Figure 6b).

3.3. Model-3

Figure 5c shows the delineated MZ map for model-3, which suggests that the central, eastern and western parts of the field have higher yield expectancy. However, a low yield pattern was observed at the edge of the field. According to Table S3, the p-value of the Kruskal–Wallis H-Test for this model was 1.2 × 10−102. Nevertheless, the pair of the zones 4–5 did not support the separability hypothesis, with p-values of 0.37 (Paired Samples T-test), 0.18 (Paired Samples T-test (log of data)) and 0.99 (Wilcoxon signed-rank test), as they are also overlapping along with the pair of the zones 1–2 (Figure 6c).

3.4. Improvement of Model Results

3.4.1. Model-1

Figure 11 shows the delineated MZ maps in 5 classes for model-1 considering agronomy and climate information. Figure 11a shows the MZ map that is based on the RS data before 25 March 2020 (see Table 2). Thus, it can be used as a fertilization recommendation map at 25 March 2020. Figure 11b shows the result of model-1 with the RS data before 21 April 2020. Additionally, this map can be utilized as a fertilizer recommendation map at 21 April 2020. As can be seen in Tables S4 and S5, all zones passed the separability hypothesis (p-value < 0.5), thus all zones are separable. The boxplots of these maps also confirmed the results of statistical tests (Figure 12). Zones 4 and 5 passed the test, though they showed partial overlap.

3.4.2. Model-3

Likewise, MZ delineation was conducted for model-3 with the hypothesis (agronomy and climate information) that was considered for improving the MZ results. Figure 13 shows the MZ maps for two dates before fertilization. As shown in Table S6, all zones were separable (p-value < 0.5). However, zones 3 and 4 did not pass the statistical test, and the p-values of statistical tests were 0.78 (Paired Samples T-test), 0.73 (Paired Samples T-test (log of data)) and 0.08 (Wilcoxon signed-rank test) (Table S7). The boxplot of both maps confirms the statistical tests (Figure 14). Similar to model-1, zones 4 and 5 were overlapping.

4. Discussion

This study followed a comparative analysis for MZ mapping that was aimed at optimizing fertilization in the context of PA. The comparative analysis has been tested by experiments based on the Sentinel-2 satellite data. The results showed only marginal consistency among the mapping outputs without any additional agronomy and climate information, which confirms the hypothesized challenges in the best model selection. Although statistical tests were considered to be appropriate validation tools, it is also suggested that the quality of the final results be eventually confirmed by a farmer and those engaged in the ongoing cultivation process.
In terms of the applied models, model-1 was partially adopted from Georgi et al. [17], and solely relied on RS data without any additional source of information. One of the advantages of this model is the utilization of NIR bands that are less noisy and more physical than band ratio vegetation indices such as NDVI. Therefore, NIR bands were used in model-3 as well. Besides, labelling of the zones from low to high was convenient compared with using a thresholding procedure for MZ delineation in other methods. On the other hand, labelling and sorting the zones in a clustering method (e.g., k-means) is a demanding task. All in all, the main disadvantage of model-1 was the lack of agronomy and climate information that can negatively impact the result.
Model-2 integrated data from multiple sources that potentially provide additional information for MZ delineation. This resulted in an enhanced accuracy of the final results. Nevertheless, one may note that such a workflow increases the costs, whereas this model is considered impractical without the incorporated sample data (soil, yield, elevation, etc.). Here, the optimum number of soil sampling data was not accessible (14 samples in 50 ha) (Table 3). Furthermore, the soil sampling data which were applied here comprised only chemical attributes and lacked soil physical attributes, such as soil compaction, texture and electric conductivity, providing further influential information on soil condition. Additionally, the certainty of utilizing RS data was confirmed in this model, as previously claimed by Martins et al. [47] and Song et al. [37]. Besides, the results of this study showed the high correlation of RS data with the yield map in correlation analysis and feature selection (See Figure 8 and Figure 9). The k-means clustering method was used in model-3 to eliminate the NIR bands’ preprocessing part (see Figure 2) and reduce the data computation of model-1. However, the labelling and sorting of the zones remain an issue in this method, as previously mentioned.
Although Georgi et al. [17] indicated that a three-year timeframe would be sufficient to depict the crop pattern for MZ delineation, the results of this study witnessed an inferior performance considering the years 2018, 2019 and 2020. Instead, focusing on one crop growing season in 2019–2020 yielded the best results. Climate change is to blame for this, as it influences crop patterns and results in year-to-year variations in crop yield. For example, in 2020, the average yield had dropped by 400 kg/ha in comparison to 2019 (see Table 3). As a result, one of the finest data sources for MZ delineation would be high-resolution remote sensing data, such as unmanned aerial vehicle (UAV) data, as it contains detailed spatial information and one UAV data point before the fertilization date can be used to create MZ and apply the generated MZ map for fertilization. Additionally, Nawar et al. [2] asserted that in terms of performance and cost, remote sensing data such as UAV or satellite data are more suitable for informing variable rate nitrogen fertilizer application than soil characteristics’ data, such as electric conductivity and soil texture.
Several studies [45,47] recently applied Cohen’s Kappa coefficient to evaluate delineated MZ maps with yield maps. However, this method was avoided due to some rationales. First, the classification of yield maps entails highly expert agronomic knowledge by which elements such as climate and agronomy information should be considered to classify the yield maps. Further, yield maps inherently comprise continuous values, subject to bias and information loss when converted and interpreted in categorical/integer data values. Last but not least, statistical tests from Georgi et al. [17] were adopted, which we think will provide an appropriate platform for validating the MZ map.
Further analyses for improvement of MZ delineation showed superior performance of model-1 and model-3 when agronomy and climate information were considered. The output MZ maps of both performed better and were similar and consistent, unlike other delineated MZ maps that lacked agronomy and climate information (see Figure 5, Figure 11 and Figure 13). This supports the importance of this information.
As an additional analysis, RS-based models were calibrated by solely considering agronomy information. Thus, the data within the phenological phases from 2018 to 2020 (2 seasons) were selected, i.e., leaving out data before seeding and after harvesting in two seasons. However, the results showed inferior performance of the RS-based models compared to the model considering both agronomy and climate information. This is parallel to suggestions presented in [4] indicating that information on crops, soil or phenological phases increases the quality of the model. Thus, incorporating both agronomy and climate information was strongly suggested in MZ delineation for PA. This issue is inevitable under new climate conditions with more frequent drought and other extreme events. As such, the requirement of considering this information is to conduct a model in a single year. Accordingly, the generated MZ map can be used as a guide for fertilization.
In a similar study, Santaga et al. [45] investigated simplified and advanced models for creating variable nitrogen fertilization maps using satellite imagery, yield maps and protein maps embedded in FMIS architecture. In contrast to the results of [45], in this study, in addition to simplified models, hybrid models were tested using non-advanced soil data, which should be considered more accessible to agricultural producers than yield or protein maps. In both studies, both simplified and hybrid or advanced models can be effectively implemented in the expanding FMIS [87].

5. Conclusions

This study compared three models based on the Sentinel-2 satellite data for management zones’ (MZ) delineation in the context of precision agriculture (PA) on an example of a winter wheat field. The remote sensing (RS)-based models 1 and 3 did not require any additional data and were thus considered less expensive than model-2. Additionally, agronomy (wheat phenological phases) and climate information (drought) were considered to improve MZ maps when applying RS-based models. The results also showed that MZ delineation is prone to uncertain results, particularly under new and rapidly changing climate conditions, if no agronomy expert knowledge and climate information were incorporated. These findings enhance the understanding of the role of new climate conditions for RS-based PA algorithms in rainfed farming with the potential agricultural drought and its impact in applying fertilizer. As drought stress is on the rise and had a significant impact on the growth of most arable crops in central Europe, it is suggested to perform low-cost RS-based techniques only for the current season. This information provides additional, reliable and current information to improve MZ delineation for optimizing fertilization in a single-year context.
In conclusion, an algorithm has been designed and its use has been evaluated for a few test cases for the integration of soil, crop and yield information, together with knowledge about agronomy and climate information. This will improve the results of MZ delineation and generate a guiding map to prescribe variable rates of fertilization before the necessary fertilizer application dates.
It is recommended for future research to use remote sensing data with high spatial resolution, such as satellite images with higher resolution and drone images, in the de-lineation of MZ maps since delineation of MZ maps with single-year data is feasible using the proposed method. It is also suggested to use cost-benefit analysis to evaluate the implemented MZ maps, e.g., in the context of variable rate fertilization. However, it was concluded that the quality of the final MZ maps should be eventually calibrated in collaboration with farmers and all these steps help farmers set up their fertilization operation to address problem areas and maximize yield.
This study has compared three models for creating MZ, focused on implementation in cloud-based farm management information systems (FMIS). The results of this work, as well as those of Santaga et al. [45], indicate that for FMIS, it is necessary to develop not only advanced models for creating MZ or variable nitrogen fertilization strategies, but also simplified and hybrid models for those users who do not have multiyear crop vegetation data or simplified soil data. Further work in this area is recommended.

Supplementary Materials

The following are available online at https://0-www-mdpi-com.brum.beds.ac.uk/article/10.3390/agriculture11111104/s1, Table S1: The results of statistical tests for model-1. The red-shaded cells show inseparability between classes (p > 0.05), Table S2: The results of statistical teste for model-2. The red-shaded colour follows the description from Table S1, Table S3: The results of statistical tests for model-3. The red-shaded colour follows the description from Table S1, Table S4: The results of statistical tests for model-1 for fertilization at 25 March 2020, Table S5: The results of statistical tests for model-1 for fertilization at 21 April 2020, Table S6: The results of statistical tests for model-3 for fertilization at 25 March 2020, Table S7: The results of statistical tests for model-3 for fertilization at 21 April 2020. The red-shaded colour follows the description from Table S1.

Author Contributions

Conceptualization, H.L., A.A.A. and M.R.; methodology, M.R., A.A.A., H.L., A.S.N. and T.W.; software, M.R.; validation, M.R., T.W., H.L. and G.N.; formal analysis, M.R., A.A.A., T.W., H.L. and G.N.; investigation, M.R., A.A.A., T.W., M.C. and H.L.; data curation, M.C., T.W., M.R., A.A.A. and G.N.; writing—original draft preparation, M.R. and H.L.; writing—review and editing, M.R., H.L., A.A.A., Y.M., T.W., A.S.N. and G.N.; visualization, M.R.; supervision, H.L., A.A.A. and Y.M. All authors have read and agreed to the published version of the manuscript.

Funding

This research has received funding from the European Union’s Horizon 2020 research and innovation programme, performed as part of the project “SmartAgriHubs—Connecting the dots to unleash the innovation potential for digital transformation of the European agri-food sector”, under grant agreement No. 818182.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The authors would like to thank KNTU for technical support. The research was conducted within the Research Laboratory “Remote Sensing for Ecology and Ecosystem Conservation (RSEEC)” of the KNTU. We would like to thank RGB Brody Poznań University of Life Sciences for providing field data. We also thank ESA for providing open access Sentinel-2 data.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Assembly, G. United Nations: Transforming Our World: The 2030 Agenda for Sustainable Development; UN: New York, NY, USA, 2015. [Google Scholar]
  2. Nawar, S.; Corstanje, R.; Halcro, G.; Mulla, D.; Mouazen, A.M. Delineation of Soil Management Zones for Variable-Rate Fertilization. Adv. Agron. 2017, 143, 175–245. [Google Scholar] [CrossRef]
  3. Foley, J.A.; Ramankutty, N.; Brauman, K.A.; Cassidy, E.S.; Gerber, J.S.; Johnston, M.; Mueller, N.D.; O’Connell, C.; Ray, D.K.; West, P.C.; et al. Solutions for a cultivated planet. Nature 2011, 478, 337–342. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  4. Sims, B.; Kienzle, J. Sustainable Agricultural Mechanization for Smallholders: What Is It and How Can We Implement It? Agriculture 2017, 7, 50. [Google Scholar] [CrossRef] [Green Version]
  5. Larkin, D.L.; Lozada, D.N.; Mason, R.E. Genomic Selection—Considerations for Successful Implementation in Wheat Breeding Programs. Agronomy 2019, 9, 479. [Google Scholar] [CrossRef] [Green Version]
  6. Gebbers, R.; Adamchuk, V.I. Precision agriculture and food security. Science 2010, 327, 828–831. [Google Scholar] [CrossRef]
  7. Lee, C.-L.; Strong, R.; Dooley, K.E. Analyzing Precision Agriculture Adoption across the Globe: A Systematic Review of Scholarship from 1999–2020. Sustainability 2021, 13, 10295. [Google Scholar] [CrossRef]
  8. Lajoie-O’Malley, A.; Bronson, K.; van der Burg, S.; Klerkx, L. The future(s) of digital agriculture and sustainable food systems: An analysis of high-level policy documents. Ecosyst. Serv. 2020, 45, 101183. [Google Scholar] [CrossRef]
  9. Pierce, F.J.; Nowak, P. Aspects of precision agriculture. Adv. Agron. 1999, 67, 1–85. [Google Scholar]
  10. Bongiovanni, R.; Lowenberg-DeBoer, J. Precision agriculture and sustainability. Precis. Agric. 2004, 5, 359–387. [Google Scholar] [CrossRef]
  11. Mueller, N.D.; Gerber, J.S.; Johnston, M.; Ray, D.K.; Ramankutty, N.; Foley, J.A. Closing yield gaps through nutrient and water management. Nature 2012, 490, 254–257. [Google Scholar] [CrossRef]
  12. Liaghat, S.; Balasundram, S.K. A review: The role of remote sensing in precision agriculture. Am. J. Agric. Biol. Sci. 2010, 5, 50–55. [Google Scholar]
  13. Wojciechowski, T.; Niedbała, G.; Czechlowski, M.; Nawrocka, J.R.; Piechnik, L.; Niemann, J. Rapeseed seeds quality classification with usage of VIS-NIR fiber optic probe and artificial neural networks. In Proceedings of the 2016 International Conference on Optoelectronics and Image Processing (ICOIP), Warsaw, Poland, 10–12 June 2016; pp. 44–48. [Google Scholar]
  14. Niazian, M.; Niedbała, G. Machine Learning for Plant Breeding and Biotechnology. Agriculture 2020, 10, 436. [Google Scholar] [CrossRef]
  15. Szwedziak, K.; Polańczyk, E.; Grzywacz, Ż.; Niedbała, G.; Wojtkiewicz, W. Neural Modeling of the Distribution of Protein, Water and Gluten in Wheat Grains during Storage. Sustainability 2020, 12, 5050. [Google Scholar] [CrossRef]
  16. Mulla, D.J. Twenty five years of remote sensing in precision agriculture: Key advances and remaining knowledge gaps. Biosyst. Eng. 2013, 114, 358–371. [Google Scholar] [CrossRef]
  17. Georgi, C.; Spengler, D.; Itzerott, S.; Kleinschmit, B. Automatic delineation algorithm for site-specific management zones based on satellite remote sensing data. Precis. Agric. 2017, 19, 684–707. [Google Scholar] [CrossRef] [Green Version]
  18. Cisternas, I.; Velásquez, I.; Caro, A.; Rodríguez, A. Systematic literature review of implementations of precision agriculture. Comput. Electron. Agric. 2020, 176, 105626. [Google Scholar]
  19. Vrindts, E.; Mouazen, A.M.; Reyniers, M.; Maertens, K.; Maleki, M.R.; Ramon, H.; De Baerdemaeker, J. Management Zones based on Correlation between Soil Compaction, Yield and Crop Data. Biosyst. Eng. 2005, 92, 419–428. [Google Scholar] [CrossRef]
  20. Fleming, K.; Heermann, D.; Westfall, D. Evaluating soil color with farmer input and apparent soil electrical conductivity for management zone delineation. Agron. J. 2004, 96, 1581–1587. [Google Scholar] [CrossRef]
  21. Van Alphen, B.; Stoorvogel, J. A Methodology to Define Management Units in Support of an Integrated, Model-Based Approach to Precision Agriculture. In Proceedings of the Fourth International Conference on Precision Agriculture, Saint Paul, MN, USA, 19–22 July 1998; pp. 1265–1278. [Google Scholar]
  22. Nolan, S.; Goddard, T.; Lohstraeter, G.; Coen, G. Assessing managements units on rolling topography. In Proceedings of the 5th International Conference on Precision Agriculture, Bloomington, MN, USA, 16–19 July, 2000; pp. 1–12. [Google Scholar]
  23. Fraisse, C.; Sudduth, K.; Kitchen, N. Delineation of site-specific management zones by unsupervised classification of topographic attributes and soil electrical conductivity. Trans. ASAE 2001, 44, 155. [Google Scholar] [CrossRef] [Green Version]
  24. Fu, Q.; Wang, Z.; Jiang, Q. Delineating soil nutrient management zones based on fuzzy clustering optimized by PSO. Math. Comput. Model. 2010, 51, 1299–1305. [Google Scholar] [CrossRef]
  25. MacMillan, R.; Pettapiece, W.; Watson, L.; Goddard, T. A landform segmentation model for precision farming. In Proceedings of the Fourth International Conference on Precision Agriculture, Saint Paul, MN, USA, 19–22 July 1998; pp. 1335–1346. [Google Scholar]
  26. Taylor, R.; Kluitenberg, G.; Schrock, M.; Zhang, N.; Schmidt, J.; Havlin, J. Using Yield Monitor Data To Determine Spati Al Crop Production Potential. Trans. ASAE 2001, 44, 1409. [Google Scholar]
  27. Diker, K.; Heermann, D.; Brodahl, M. Frequency analysis of yield for delineating yield response zones. Precis. Agric. 2004, 5, 435–444. [Google Scholar] [CrossRef]
  28. Milne, A.E.; Webster, R.; Ginsburg, D.; Kindred, D. Spatial multivariate classification of an arable field into compact management zones based on past crop yields. Comput. Electron. Agric. 2012, 80, 17–30. [Google Scholar] [CrossRef]
  29. Niedbała, G.; Kozlowski, J. Application of artificial neural networks for multi-criteria yield prediction of winter wheat. J. Agric. Sci. Technol. 2019, 21, 51–61. [Google Scholar]
  30. Niedbała, G.; Nowakowski, K.; Rudowicz-Nawrocka, J.; Piekutowska, M.; Weres, J.; Tomczak, R.J.; Tyksiński, T.; Álvarez Pinto, A. Multicriteria prediction and simulation of winter wheat yield using extended qualitative and quantitative data based on artificial neural networks. Appl. Sci. 2019, 9, 2773. [Google Scholar] [CrossRef] [Green Version]
  31. Hara, P.; Piekutowska, M.; Niedbała, G. Selection of Independent Variables for Crop Yield Prediction Using Artificial Neural Network Models with Remote Sensing Data. Land 2021, 10, 609. [Google Scholar] [CrossRef]
  32. Kitchen, N.R.; Sudduth, K.A.; Myers, D.B.; Drummond, S.T.; Hong, S.Y. Delineating productivity zones on claypan soil fields using apparent soil electrical conductivity. Comput. Electron. Agric. 2005, 46, 285–308. [Google Scholar] [CrossRef]
  33. Cambouris, A.; Nolin, M.; Zebarth, B.; Laverdière, M. Soil management zones delineated by electrical conductivity to characterize spatial and temporal variations in potato yield and in soil properties. Am. J. Potato Res. 2006, 83, 381–395. [Google Scholar]
  34. Ahn, C.W.; Baumgardner, M.; Biehl, L. Delineation of soil variability using geostatistics and fuzzy clustering analyses of hyperspectral data. Soil Sci. Soc. Am. J. 1999, 63, 142–150. [Google Scholar]
  35. Vizzari, M.; Santaga, F.; Benincasa, P. Sentinel 2-Based Nitrogen VRT Fertilization in Wheat: Comparison between Traditional and Simple Precision Practices. Agronomy 2019, 9, 278. [Google Scholar] [CrossRef] [Green Version]
  36. Fridgen, J.J.; Kitchen, N.R.; Sudduth, K.A.; Drummond, S.T.; Wiebold, W.J.; Fraisse, C.W. Management Zone Analyst (MZA) Software for Subfield Management Zone Delineation. Agron. J. 2004, 96, 100–108. [Google Scholar]
  37. Song, X.; Wang, J.; Huang, W.; Liu, L.; Yan, G.; Pu, R. The delineation of agricultural management zones with high resolution remotely sensed data. Precis. Agric. 2009, 10, 471–487. [Google Scholar] [CrossRef]
  38. De Benedetto, D.; Castrignanò, A.; Rinaldi, M.; Ruggieri, S.; Santoro, F.; Figorito, B.; Gualano, S.; Diacono, M.; Tamborrino, R. An approach for delineating homogeneous zones by using multi-sensor data. Geoderma 2013, 199, 117–127. [Google Scholar] [CrossRef]
  39. Yao, R.-J.; Yang, J.-S.; Zhang, T.-J.; Gao, P.; Wang, X.-P.; Hong, L.-Z.; Wang, M.-W. Determination of site-specific management zones using soil physico-chemical properties and crop yields in coastal reclaimed farmland. Geoderma 2014, 232–234, 381–393. [Google Scholar] [CrossRef]
  40. Farid, H.U.; Bakhsh, A.; Ahmad, N.; Ahmad, A.; Mahmood-Khan, Z. Delineating site-specific management zones for precision agriculture. J. Agric. Sci. 2015, 154, 273–286. [Google Scholar] [CrossRef]
  41. Peralta, N.R.; Costa, J.L.; Balzarini, M.; Castro Franco, M.; Córdoba, M.; Bullock, D. Delineation of management zones to improve nitrogen management of wheat. Comput. Electron. Agric. 2015, 110, 103–113. [Google Scholar] [CrossRef]
  42. Basso, B.; Dumont, B.; Cammarano, D.; Pezzuolo, A.; Marinello, F.; Sartori, L. Environmental and economic benefits of variable rate nitrogen fertilization in a nitrate vulnerable zone. Sci. Total Environ. 2016, 545–546, 227–235. [Google Scholar] [CrossRef] [Green Version]
  43. Gavioli, A.; de Souza, E.G.; Bazzi, C.L.; Guedes, L.P.C.; Schenatto, K. Optimization of management zone delineation by using spatial principal components. Comput. Electron. Agric. 2016, 127, 302–310. [Google Scholar] [CrossRef]
  44. Oldoni, H.; Silva Terra, V.S.; Timm, L.C.; Júnior, C.R.; Monteiro, A.B. Delineation of management zones in a peach orchard using multivariate and geostatistical analyses. Soil Tillage Res. 2019, 191, 1–10. [Google Scholar] [CrossRef]
  45. Gavioli, A.; de Souza, E.G.; Bazzi, C.L.; Schenatto, K.; Betzek, N.M. Identification of management zones in precision agriculture: An evaluation of alternative cluster analysis methods. Biosyst. Eng. 2019, 181, 86–102. [Google Scholar] [CrossRef]
  46. Moharana, P.C.; Jena, R.K.; Pradhan, U.K.; Nogiya, M.; Tailor, B.L.; Singh, R.S.; Singh, S.K. Geostatistical and fuzzy clustering approach for delineation of site-specific management zones and yield-limiting factors in irrigated hot arid environment of India. Precis. Agric. 2019, 21, 426–448. [Google Scholar] [CrossRef]
  47. Nogueira Martins, R.; Magalhães Valente, D.S.; Fim Rosas, J.T.; Souza Santos, F.; Lima Dos Santos, F.F.; Nascimento, M.; Campana Nascimento, A.C. Site-specific Nutrient Management Zones in Soybean Field Using Multivariate Analysis: An Approach Based on Variable Rate Fertilization. Commun. Soil Sci. Plant. Anal. 2020, 51, 687–700. [Google Scholar] [CrossRef]
  48. Gotway, C.; Ferguson, R.; Hergert, G. The Effects of Mapping and Scale on Variable-Rate Fertilizer Recommendations for Corn. In Proceedings of the Third International Conference on Precision Agriculture, Minneapolis, MN, USA, 23–26 June 1996; pp. 321–330. [Google Scholar]
  49. Fleming, K.; Westfall, D.; Bausch, W. Evaluating management zone technology and grid soil sampling for variable rate nitrogen application. In Proceedings of the 5th International Conference on Precision Agriculture, Bloomington, MN, USA, 16–19 July 2000; pp. 16–19. [Google Scholar]
  50. Koch, B.; Khosla, R. The Role of Precision Agriculture in Cropping Systems. J. Crop. Prod. 2003, 9, 361–381. [Google Scholar] [CrossRef]
  51. Cohen, S.; Cohen, Y.; Alchanatis, V.; Levi, O. Combining spectral and spatial information from aerial hyperspectral images for delineating homogenous management zones. Biosyst. Eng. 2013, 114, 435–443. [Google Scholar] [CrossRef]
  52. Basnyat, P.; McConkey, B.G.; Selles, F.; Meinert, L.B. Effectiveness of using vegetation index to delineate zones of different soil and crop grain production characteristics. Can. J. Soil Sci. 2005, 85, 319–328. [Google Scholar] [CrossRef]
  53. Oza, S.R.; Panigrahy, S.; Parihar, J.S. Concurrent use of active and passive microwave remote sensing data for monitoring of rice crop. Int. J. Appl. Earth Obs. Geoinf. 2008, 10, 296–304. [Google Scholar] [CrossRef]
  54. Segarra, J.; Buchaillot, M.L.; Araus, J.L.; Kefauver, S.C. Remote Sensing for Precision Agriculture: Sentinel-2 Improved Features and Applications. Agronomy 2020, 10, 641. [Google Scholar] [CrossRef]
  55. European Space Agency (ESA). MultiSpectral Instrument (MSI) Overview. Available online: https://sentinels.copernicus.eu/web/sentinel/technical-guides/sentinel-2-msi/msi-instrument (accessed on 20 October 2021).
  56. European Space Agency (ESA). Sentinel-2 (MSI) Overview. Available online: https://sentinel.esa.int/web/sentinel/missions/sentinel-2 (accessed on 20 October 2021).
  57. European Space Agency (ESA). Data Products (MSI) Overview. Available online: https://sentinels.copernicus.eu/web/sentinel/missions/sentinel-2/data-products (accessed on 20 October 2021).
  58. Rouse, J.W.; Haas, R.H.; Schell, J.A.; Deering, D.W. Monitoring vegetation systems in the Great Plains with ERTS. NASA Spec. Publ. 1974, 351, 309. [Google Scholar]
  59. Weiss, M.; Baret, F. S2ToolBox Level 2 Products: LAI, FAPAR, FCOVER, Version 1.1. Available online: https://step.esa.int/docs/extra/ATBD_S2ToolBox_L2B_V1.1.pdf (accessed on 20 October 2021).
  60. Song, W.; Mu, X.; Ruan, G.; Gao, Z.; Li, L.; Yan, G. Estimating fractional vegetation cover and the vegetation index of bare soil and highly dense vegetation with a physically based method. Int. J. Appl. Earth Obs. Geoinf. 2017, 58, 168–176. [Google Scholar]
  61. Kosior, K. Łańcuch wartości dużych zbiorów danych (big data) w rolnictwie–problemy i wyzwania regulacyjne. Studia BAS 2020, 3, 101–125. [Google Scholar] [CrossRef]
  62. Santaga, F.S.; Benincasa, P.; Toscano, P.; Antognelli, S.; Ranieri, E.; Vizzari, M. Simplified and Advanced Sentinel-2-Based Precision Nitrogen Management of Wheat. Agronomy 2021, 11, 1156. [Google Scholar]
  63. Majchrzak, L.; Sawinska, Z.; Natywa, M.; Skrzypczak, G.; Głowicka-Wołoszyn, R. Impact of different tillage systems on soil dehydrogenase activity and spring wheat infection. J. Agric. Sci. Technol. 2016, 18, 1871–1881. [Google Scholar]
  64. Tryjanowski, P.; Sparks, T.H.; Blecharczyk, A.; Małecka-Jankowiak, I.; Switek, S.; Sawinska, Z. Changing Phenology of Potato and of the Treatment for its Major Pest (Colorado Potato Beetle)—A Long-term Analysis. Am. J. Potato Res. 2017, 95, 26–32. [Google Scholar] [CrossRef] [Green Version]
  65. Egnér, H.; Riehm, H.; Domingo, W. Untersuchungen über die chemische Bodenanalyse als Grundlage für die Beurteilung des Nährstoffzustandes der Böden. II. Chem. Extra Ktionsmethoden Zur Phosphor-Und Kaliumbestimmung K. Lantbr. Ann. 1960, 26, 199–215. [Google Scholar]
  66. Schachtschabel, P.V. Das pflanzenverfügbare Magnesium des Boden und seine Bestimmung. Z. Für Pflanz. Düngung Bodenkd. 1954, 67, 9–23. [Google Scholar] [CrossRef]
  67. Tiurin, I. К mietodikie analiza dla srawnitielnogo izuczenia sostawa poczwiennogo gumusa. Tr. Poczw. Inst. AN SSSR 1951, 38, 1–250. [Google Scholar]
  68. Czechlowski, M.; Wojciechowski, T. The utilization of information about local variable environmental conditions to predict the quality of wheat grain during the harvest. J. Res. Appl. Agric. Eng. 2013, 58, 31–34. [Google Scholar]
  69. Czechlowski, M.; Wojciechowski, T. System mechatroniczny do selektywnego zbioru ziarna zbóż. Zeszyty Naukowe Instytutu Pojazdów 2013, 4, 95. [Google Scholar]
  70. Czechlowski, M.; Wojciechowski, T.; Adamski, M.; Niedbała, G.; Piekutowska, M. Application of ASG-EUPOS high precision positioning system for cereal harvester monitoring. J. Res. Appl. Agric. Eng. 2018, 63, 44–50. [Google Scholar]
  71. European Space Agency (ESA). Copernicus Open Access Hub. Available online: https://scihub.copernicus.eu/ (accessed on 20 October 2021).
  72. European Space Agency (ESA). STEP-Scientific Toolbox Exploitation Platform Ver 7.0. Available online: https://step.esa.int/main/ (accessed on 20 October 2021).
  73. Liang, S. Quantitative Remote Sensing of Land Surfaces; John Wiley & Sons Inc.: Hoboken, NJ, USA, 2005; Volume 30. [Google Scholar]
  74. Lillesand, T.; Kiefer, R.W.; Chipman, J. Remote Sensing and Image Interpretation; John Wiley & Sons Inc.: Hoboken, NJ, USA, 2015. [Google Scholar]
  75. QGIS, Version 3.10.8. Available online: https://www.qgis.org/en/site/ (accessed on 8 February 2021).
  76. Webster, R.; Oliver, M.A. Geostatistics for Environmental Scientists; John Wiley & Sons Inc.: Hoboken, NJ, USA, 2007. [Google Scholar]
  77. ArcGIS, Version 10.7. Available online: https://desktop.arcgis.com/en/arcmap/ (accessed on 8 February 2021).
  78. Management Zone Analyst (MZA), Version 1.0. Available online: https://www.ars.usda.gov/research/software/download/?softwareid=24&modecode=50-70-10-00 (accessed on 8 February 2021).
  79. Odeh, I.; McBratney, A.; Chittleborough, D. Soil pattern recognition with fuzzy-c-means: Application to classification and soil-landform interrelationships. Soil Sci. Soc. Am. J. 1992, 56, 505–516. [Google Scholar] [CrossRef]
  80. Boydell, B.; McBratney, A. Identifying potential within-field management zones from cotton-yield estimates. Precis. Agric. 2002, 3, 9–23. [Google Scholar]
  81. Bezdek, J.C. Pattern Recognition with Fuzzy Objective Function Algorithms; Springer: Boston, MA, USA, 1981; pp. 43–93. [Google Scholar]
  82. Olofsson, P.; Foody, G.M.; Herold, M.; Stehman, S.V.; Woodcock, C.E.; Wulder, M.A. Good practices for estimating area and assessing accuracy of land change. Remote. Sens. Environ. 2014, 148, 42–57. [Google Scholar] [CrossRef]
  83. Pińskwar, I.; Choryński, A.; Kundzewicz, Z.W. Severe Drought in the Spring of 2020 in Poland—More of the Same? Agronomy 2020, 10, 1646. [Google Scholar] [CrossRef]
  84. Warrick, A. Spatial variability of soil physical properties in the field. Appl. Soil Phys. 1980, 319–344. [Google Scholar] [CrossRef]
  85. Cambardella, C.A.; Moorman, T.B.; Novak, J.; Parkin, T.; Karlen, D.; Turco, R.; Konopka, A. Field-scale variability of soil properties in central Iowa soils. Soil Sci. Soc. Am. J. 1994, 58, 1501–1511. [Google Scholar] [CrossRef]
  86. Reza, S.K.; Nayak, D.C.; Mukhopadhyay, S.; Chattopadhyay, T.; Singh, S.K. Characterizing spatial variability of soil properties in alluvial soils of India using geostatistics and geographical information system. Arch. Agron. Soil Sci. 2017, 63, 1489–1498. [Google Scholar] [CrossRef]
  87. Giua, C.; Materia, V.C.; Camanzi, L. Management information system adoption at the farm level: Evidence from the literature. Br. Food J. 2020, 123, 884–909. [Google Scholar] [CrossRef]
Agriculture 11 01104 g001 550
Figure 1. The geographic location of the study area in Brody, Wielkopolska province, Poland.
Figure 1. The geographic location of the study area in Brody, Wielkopolska province, Poland.
Agriculture 11 01104 g001
Agriculture 11 01104 g002 550
Figure 2. The workflow of model-1 (RS- and threshold-based clustering).
Figure 2. The workflow of model-1 (RS- and threshold-based clustering).
Agriculture 11 01104 g002
Agriculture 11 01104 g003 550
Figure 3. The workflow of model-2 (hybrid-based, unsupervised clustering).
Figure 3. The workflow of model-2 (hybrid-based, unsupervised clustering).
Agriculture 11 01104 g003
Agriculture 11 01104 g004 550
Figure 4. The workflow of model-3 (RS-based, unsupervised clustering).
Figure 4. The workflow of model-3 (RS-based, unsupervised clustering).
Agriculture 11 01104 g004
Agriculture 11 01104 g005 550
Figure 5. Delineated MZ maps of model-1 (a), model-2 (b) and model-3 (c).
Figure 5. Delineated MZ maps of model-1 (a), model-2 (b) and model-3 (c).
Agriculture 11 01104 g005
Agriculture 11 01104 g006 550
Figure 6. Boxplots of stratified sampling for model-1 (a), model-2 (b) and model-3 (c). The white circles out of the whiskers show outliers.
Figure 6. Boxplots of stratified sampling for model-1 (a), model-2 (b) and model-3 (c). The white circles out of the whiskers show outliers.
Agriculture 11 01104 g006
Agriculture 11 01104 g007 550
Figure 7. Spatial distribution maps of soil attributes including pHKCl (a), P2O5 (b), K2O (c), Mg (d) and OM (e), the spatial distribution maps of RS data including NDVI (f), LAI (g), FCV (h) and FAPAR (i), the spatial distribution maps of yield for 2019 (j) and 2020 (k) and the spatial distribution maps of geomorphology data including DEM (l) and Slope (m).
Figure 7. Spatial distribution maps of soil attributes including pHKCl (a), P2O5 (b), K2O (c), Mg (d) and OM (e), the spatial distribution maps of RS data including NDVI (f), LAI (g), FCV (h) and FAPAR (i), the spatial distribution maps of yield for 2019 (j) and 2020 (k) and the spatial distribution maps of geomorphology data including DEM (l) and Slope (m).
Agriculture 11 01104 g007
Agriculture 11 01104 g008 550
Figure 8. Spearman’s correlation matrix amongst soil attributes, RS data, yield data and geomorphology data.
Figure 8. Spearman’s correlation matrix amongst soil attributes, RS data, yield data and geomorphology data.
Agriculture 11 01104 g008
Agriculture 11 01104 g009 550
Figure 9. Result of RF feature importance.
Figure 9. Result of RF feature importance.
Agriculture 11 01104 g009
Agriculture 11 01104 g010 550
Figure 10. Fuzzy performance index (FPI) and normalized classification entropy (NCE) calculated to determine the optimum number of clusters for the study area.
Figure 10. Fuzzy performance index (FPI) and normalized classification entropy (NCE) calculated to determine the optimum number of clusters for the study area.
Agriculture 11 01104 g010
Agriculture 11 01104 g011 550
Figure 11. Delineated MZ maps of model-1 for fertilization at 25 March 2020 (a) and 21 April 2020 (b).
Figure 11. Delineated MZ maps of model-1 for fertilization at 25 March 2020 (a) and 21 April 2020 (b).
Agriculture 11 01104 g011
Agriculture 11 01104 g012 550
Figure 12. Boxplots of stratified sampling for model-1 for fertilizations at 25 March 2020 (a) and 21 April 2020 (b). The white circles out of the whiskers show outliers.
Figure 12. Boxplots of stratified sampling for model-1 for fertilizations at 25 March 2020 (a) and 21 April 2020 (b). The white circles out of the whiskers show outliers.
Agriculture 11 01104 g012
Agriculture 11 01104 g013 550
Figure 13. Delineated MZ maps of model-3 for fertilization at 25 March 2020 (a) and 21 April 2020 (b).
Figure 13. Delineated MZ maps of model-3 for fertilization at 25 March 2020 (a) and 21 April 2020 (b).
Agriculture 11 01104 g013
Agriculture 11 01104 g014 550
Figure 14. Boxplots of stratified sampling for model-3 for fertilization at 25 March 2020 (a) and 21 April 2020 (b). The white circles out of the whiskers show outliers.
Figure 14. Boxplots of stratified sampling for model-3 for fertilization at 25 March 2020 (a) and 21 April 2020 (b). The white circles out of the whiskers show outliers.
Agriculture 11 01104 g014
Table
Table 1. Monthly distribution of 119 Sentinel-2 L2A scenes available from 1 January 2018 to 1 July 2020.
Table 1. Monthly distribution of 119 Sentinel-2 L2A scenes available from 1 January 2018 to 1 July 2020.
MonthJanFebMarAprMayJuneJulyAugSepOctNovDec
Year
2018121582675821
2019042658335524
2020334743- *- *- *- *- *- *
* Not analysed.
Table
Table 2. Acquisition dates of final selected Sentinel-2 data.
Table 2. Acquisition dates of final selected Sentinel-2 data.
DataDataData
25 February 20188 February 201920 February 2020
6 April 201818 February 201911 March 2020
9 April 201825 February 20195 April 2020
31 May 201828 February 20198 April 2020
3 June 201818 June 201922 June 2020
8 June 201820 June 2019
3 July 201825 June 2019
31 October 201824 August 2019
7 November 201827 August 2019
5 December 2018
Total in 2018: 10Total in 2019: 9Total in 2020: 5
Table
Table 3. Descriptive statistics of DEM, soil attributes and wheat yield.
Table 3. Descriptive statistics of DEM, soil attributes and wheat yield.
AttributesnMinMaxMeanSDSECVSkewnessKurtosis
DEM 2020 (m)67,15892.42100.1296.212.360.010.020.06−1.45
OM 2011 (%)521.142.641.660.440.060.260.92−0.33
pHKCl 2020146.007.16.440.320.080.050.44−0.45
P2O 2020 (mg 100 g soil−1)1417.2036.625.715.921.580.220.20−1.07
K2O 2020 (mg 100 g soil−1)1423.0034.026.863.550.950.130.64−0.80
Mg 2020 (mg 100 g soil−1)148.5012.49.891.410.380.140.70−1.04
Yield 2019 (t ha−1)96130.3915.187.191.560.020.22−1.255.01
Yield 2020 (t ha−1)85200.3613.126.821.580.020.23−1.303.11
Table
Table 4. Semi-variogram parameters of DEM, soil attributes and wheat yields.
Table 4. Semi-variogram parameters of DEM, soil attributes and wheat yields.
VariablesModelNugget
(C0)		Partial Sill
(C1)		Sill
(C0 + C1)		Nugget/Sill
C0/(C0 + C1)		Range (m)RMSE
DEMExponential00.00050.000501.40420.0209
OMJ-Bessel0.02660.24430.27090.09821247.50.1665
pHKClJ-Bessel0.02990.07780.10770.2776925.030.2391
P2O5Hole Effect12.88328.33141.2140.3126915.834.2782
K2O Gaussian6.456016.51622.9720.28101147.22.8914
MgGaussian1.01102.73383.74480.29701147.21.0427
Yield 2019Exponential1.61690.73822.35510.6865490.051.3108
Yield 2020Exponential1.22262.22083.44340.35501301.71.0925
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Rokhafrouz, M.; Latifi, H.; Abkar, A.A.; Wojciechowski, T.; Czechlowski, M.; Naieni, A.S.; Maghsoudi, Y.; Niedbała, G. Simplified and Hybrid Remote Sensing-Based Delineation of Management Zones for Nitrogen Variable Rate Application in Wheat. Agriculture 2021, 11, 1104. https://0-doi-org.brum.beds.ac.uk/10.3390/agriculture11111104

AMA Style

Rokhafrouz M, Latifi H, Abkar AA, Wojciechowski T, Czechlowski M, Naieni AS, Maghsoudi Y, Niedbała G. Simplified and Hybrid Remote Sensing-Based Delineation of Management Zones for Nitrogen Variable Rate Application in Wheat. Agriculture. 2021; 11(11):1104. https://0-doi-org.brum.beds.ac.uk/10.3390/agriculture11111104

Chicago/Turabian Style

Rokhafrouz, Mohammad, Hooman Latifi, Ali A. Abkar, Tomasz Wojciechowski, Mirosław Czechlowski, Ali Sadeghi Naieni, Yasser Maghsoudi, and Gniewko Niedbała. 2021. "Simplified and Hybrid Remote Sensing-Based Delineation of Management Zones for Nitrogen Variable Rate Application in Wheat" Agriculture 11, no. 11: 1104. https://0-doi-org.brum.beds.ac.uk/10.3390/agriculture11111104

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop