Assessing Soil and Crop Characteristics at Sub-Field Level Using Unmanned Aerial System and Geospatial Analysis

: Practicing agriculture is a multiparametric and for this reason demanding task. It involves the management of many factors and thorough strategic planning in a highly variable and uncertain environment. Crop production is a function of agricultural practices as applied in natural resources, such as soil and plants. support system, it can be used by policy makers for adapting circular economy strategies in crop production.


Introduction
A physical environment has a dynamic behavior, as it is an open system and involves variable driving forces that interact with each other in a spatial and temporal scale [1]. An anthropogenic environment is another dynamic system, which intervenes to the functions of the physical environment and highly influences its elements. A typical example of recovery, and under special treatment, it can lead to diminished waste load [23]. Regarding waste management, a comprehensive framework for conceptualizing and implementing a strategy is analyzed in [24], while the evaluation of environmental performance with the use of key indicators is realized in [25], which may also provide insights and guide strategic planning in crop management. Soil-aquifer pollution and agricultural waste strategies receive great attention in literature [26][27][28] while site-specific crop management may positively contribute to applying those managerial plans in an effective and sustainable way [29].
Remote sensing technology is a very common means, not only for monitoring crop variability [30][31][32], but also for environmental applications in general [33]. The evolution of small, unmanned aerial vehicles (UAVs) integrated with a variety of payload sensing systems, which are available in low-cost and mass production hardware solutions, has led to their wide use in agricultural research experimentation [34]. Indeed, many aerial datasets have been acquired via UAVs, and analyzed, correlating spectral signatures with parameters of field conditions [35]. The spatial nature of these datasets, along with the point reference of soil and crop field measurements, make the use of geospatial methodologies apparent, as far as their analysis is concerned [36]. Geographic information systems (GIS) offer an integrated suite of software and hardware tools for managing spatiotemporal data and mapping in-field variability [37,38]. Many research works have combined the aforementioned methods, mostly for evaluating soil conditions [39][40][41], while monitoring in situ crop nutrient variations is not so commonly met in scientific literature. Several studies have also used geospatial analysis for separating a field into management zones, concluding that this practice can well manage sub-field soil and crop variability [42][43][44].
The main purpose of conducting the present study was to investigate possible interrelations between spectral and field data under the concept of site-specific crop management. Towards this objective, a commercial low-cost UAV was used to capture orthophotographs in near-infrared wavelength and classic field sampling and analyses procedures were incorporated to determine soil and crop parameters in cotton crops. Consequently, statistical and geostatistical methods combined with the fuzzy algorithm application, undertook the correlation and the assessment of the outputs. A further scope was to check the feasibility of using a UAV in order to delineate in-field management zones based on soil, crop, and reflectance data.

Study Area
The study area was located in Chaeronea (38 • 29 44 N, 22 • 50 41 E), central Greece. Several visits to the area were conducted in order to visually inspect the cultivated land, meet local farmers, and further select the appropriate fields for the experimentation. Criteria for the selection of the fields were in-field soil, topographical variations, and favorable accessibility conditions for fieldwork. The final selection concerned two fields sown with cotton. In both pilot fields, the cotton variety sown was Celia, the distance between the rows was 95 cm, and the sowing density was 23 to 25 plants per meter. The area of the first field (Field 1) was 2.3 ha and it was completely flat, while the second one (Field 2) was a 1.9 ha sloped field ( Figure 1). The experimental plan involved one growing period.

Data Acquisition
Data acquisition concerned a series of field operations that took place throughout an entire cultivation period of the two fields and concerned soil sampling, leaf sampling, and aerial image capturing.

Soil Samples
Preliminarily to soil sampling, a grid plan was compiled for each field based on two aerial photos of the fields, with the use of the GIS software package ArcGIS (ESRI Inc., California, USA). Grid pixel size was 25 m and, hence, the number of sampling positions was 45 in Field 1 and 50 in Field 2 ( Figure 2 and 3). In order to spatially cover as representatively as possible, the area of the two fields, several in-field sampling sites exceeded the grid boundaries. Soil sampling took place in April with the use of a standard soil auger. The depth of the soil samples was 0-30 cm. Soil physicochemical analyses included the measurement of soil texture (clay content), soil acidity (pH), calcium carbonate percentage (CaCO3), organic matter (OM), total nitrogen, and electrical conductivity (EC). The descriptive statistics of the aforementioned soil parameters per field are presented in Table 1

Data Acquisition
Data acquisition concerned a series of field operations that took place throughout an entire cultivation period of the two fields and concerned soil sampling, leaf sampling, and aerial image capturing.

Soil Samples
Preliminarily to soil sampling, a grid plan was compiled for each field based on two aerial photos of the fields, with the use of the GIS software package ArcGIS (ESRI Inc., California, USA). Grid pixel size was 25 m and, hence, the number of sampling positions was 45 in Field 1 and 50 in Field 2 ( Figures 2 and 3). In order to spatially cover as representatively as possible, the area of the two fields, several in-field sampling sites exceeded the grid boundaries. Soil sampling took place in April with the use of a standard soil auger. The depth of the soil samples was 0-30 cm. Soil physicochemical analyses included the measurement of soil texture (clay content), soil acidity (pH), calcium carbonate percentage (CaCO 3 ), organic matter (OM), total nitrogen, and electrical conductivity (EC). The descriptive statistics of the aforementioned soil parameters per field are presented in Tables 1 and 2.

Leaf Samples
Leaf sampling in the study fields took place in July. Both cotton fields were at the full bloom stage, and each sample contained approximately 30 late mature leaves cut from different plants covering an area of one square meter around the sampling site. The number and the actual position of the sampling sites were based on the grid plan originally drawn for soil sampling, and regarded a subset of the respective soil samples. In each field, leaf samples were collected from 30 sites (Figures 4 and 5). The leaves after cutting were put in paper bags and transferred to the laboratory for further analyses. The leaf samples were analyzed for total nitrogen, potassium, iron, copper, and zinc. The descriptive statistics of the macro and micronutrients measured in the cotton leaves in the two fields are presented in Tables 3 and 4.  Leaf sampling in the study fields took place in July. Both cotton fields were at the full bloom stage, and each sample contained approximately 30 late mature leaves cut from different plants covering an area of one square meter around the sampling site. The number and the actual position of the sampling sites were based on the grid plan originally drawn for soil sampling, and regarded a subset of the respective soil samples. In each field, leaf samples were collected from 30 sites (Figure 4 and 5). The leaves after cutting were put in paper bags and transferred to the laboratory for further analyses. The leaf samples were analyzed for total nitrogen, potassium, iron, copper, and zinc. The descriptive statistics of the macro and micronutrients measured in the cotton leaves in the two fields are presented in Table 3; Table 4.    Aerial data were captured with a UAS. The UAS comprised of a low-cost and lowweight quadcopter (Phantom, Da-Jiang Innovations (DJI) Science and Technology Co., Ltd., Shenzhen, China) and a 12-Megapixel (4000 × 3000 pixels) Hero digital camera (GoPro Inc., California, USA), especially modified by Back-Bone Gear Inc., Ottawa Canada ( Figure 6).

Aerial Data
Aerial data were captured with a UAS. The UAS comprised of a low-cost and lowweight quadcopter (Phantom, Da-Jiang Innovations (DJI) Science and Technology Co., Ltd., Shenzhen, China) and a 12-Megapixel (4000 × 3000 pixels) Hero digital camera (GoPro Inc., San Mateo, CA, USA), especially modified by Back-Bone Gear Inc., Ottawa Canada ( Figure 6). The camera's stock lens was substituted by a flat lens (focal length 4.4 mm, relative aperture f/2.8, sensor size 1/2.3 inch and 71° angle of view). This modification was made in order to eliminate the distortion produced by fish-eye effect of the wide-angle lens [45]. Moreover, the modified camera was configured with removable filter positioning. Particularly, two detachable filters (IR-cut [400-700 nm] and near-infrared long-pass [715-1100 nm]) were used to capture aerial photographs in visible and near-infrared wavelengths, respectively.
The digital camera was attached to the quadcopter on a gimbal for compensating the quadcopter movement during flights and guaranteeing nadir image acquisition [46]. The camera was triggered before each flight mission to capture one image per 2 s. The quadcopter was equipped with the standard DJI flight control system, consisting of the main controller (MC), the inertial measurement unit (IMU), the global positioning system (GPS) receiver, and the compass. The quadcopter was remotely piloted through a 2.4 GHz radio control transmitter. Two flight missions per pilot field were scheduled and carried out. The first one was in April and the second in July, as close as possible to the dates of soil and plant sampling, respectively. The flights were carried out in stable cloud free ambient light conditions in the morning, at about 11:00 a.m. local time, at 250 m above ground level (AGL). The data collected by the UAS were two pairs of aerial images per pilot field and per flight date. Each pair consisted of one color composite (red-green-blue) image and one image in near-infrared wavelength. The selection of taking one image per field was made because photo stitching or orthomosaic generation would have introduced errors and artifacts in aerial data. The parameters used in all further analyses were the digital numbers (DN) of the aerial photos as captured in the near-infrared area of spectrum [47,48].

Data Management and Analysis
The soil parameters used in the analysis for Field 1 were those presented in Table 1, while for the analysis of Field 2, slope parameter (%) was included. Slope was calculated using an existing digital surface model (DSM) of the area. Respectively, the list of leaf The camera's stock lens was substituted by a flat lens (focal length 4.4 mm, relative aperture f/2.8, sensor size 1/2.3 inch and 71 • angle of view). This modification was made in order to eliminate the distortion produced by fish-eye effect of the wide-angle lens [45]. Moreover, the modified camera was configured with removable filter positioning. Particularly, two detachable filters (IR-cut [400-700 nm] and near-infrared long-pass [715-1100 nm]) were used to capture aerial photographs in visible and near-infrared wavelengths, respectively.
The digital camera was attached to the quadcopter on a gimbal for compensating the quadcopter movement during flights and guaranteeing nadir image acquisition [46]. The camera was triggered before each flight mission to capture one image per 2 s. The quadcopter was equipped with the standard DJI flight control system, consisting of the main controller (MC), the inertial measurement unit (IMU), the global positioning system (GPS) receiver, and the compass. The quadcopter was remotely piloted through a 2.4 GHz radio control transmitter. Two flight missions per pilot field were scheduled and carried out. The first one was in April and the second in July, as close as possible to the dates of soil and plant sampling, respectively. The flights were carried out in stable cloud free ambient light conditions in the morning, at about 11:00 a.m. local time, at 250 m above ground level (AGL). The data collected by the UAS were two pairs of aerial images per pilot field and per flight date. Each pair consisted of one color composite (red-green-blue) image and one image in near-infrared wavelength. The selection of taking one image per field was made because photo stitching or orthomosaic generation would have introduced errors and artifacts in aerial data. The parameters used in all further analyses were the digital numbers (DN) of the aerial photos as captured in the near-infrared area of spectrum [47,48].

Data Management and Analysis
The soil parameters used in the analysis for Field 1 were those presented in Table 1, while for the analysis of Field 2, slope parameter (%) was included. Slope was calculated using an existing digital surface model (DSM) of the area. Respectively, the list of leaf parameters used is shown in Table 3. All data were georeferenced and stored in a geodatabase (ArcGIS, ESRI Inc., Redlands, CA, USA) for further processing. Data management concerned the harmonization of formats, meaning that all data were transformed into Sustainability 2021, 13, 2855 9 of 24 raster files. Point referenced data were converted into continuous raster surfaces with the use of ordinary kriging geostatistical method of interpolation [49]. Aerial imagery yet being in raster format needed no further transformation. Nevertheless, in order to correctly (and accurately) apply overlay spatial methods to all datasets, common pixel size (0.1 m) and grid size, starting point, and orientation were met. The raster files, after being overlaid, were extracted to a single attribute table per field, in which every row represented a pixel and every column contained the values of the parameters involved in each analysis.
The succeeded steps of the analysis were accomplished through a combination of principal component analysis (PCA) and of Fuzzy c-means algorithm application. PCA is a multivariate statistical technique that aims to reduce the dimensions of a multiparameter system, by recognizing linear combinations between its initial variables. It aims to not only reduce the number of initial data, but also to analyze and explore interrelations between them [50]. More specifically, PCA was executed using SPSS (IBM Corp., Armonk, NY, USA) in the attribute data of the two fields, in order to investigate possible relationships or groupings between those parameters. In the applied PCA, Varimax with Kaiser normalization method of rotation was selected, and Kaiser's rule was decided for preserving the principal components [51,52]. PCA scores, accounting for every pixel of the combined raster file, underwent fuzzy classification using the Fuzzy c-means algorithm as applied in MATLAB software (MathWorks Inc., Natick, MA, USA). Fuzzy c-means classifies each participating value into a selected number of clusters with variable membership values ranging between 0 and 1. It is an unsupervised iterative algorithm that finds the cluster centroids of groups for which the in-cluster variability is minimized, whereas the between-cluster variability is maximized [53]. A flowchart of the methodological steps followed in the analysis is presented in Figure 7. The analysis shown in Figure 7 was followed twice for each field. Regarding Field 1, the first one was using soil and aerial data (April flight) and the second using crop and aerial data (July flight). Similarly, for Field 2 the analysis was run firstly with soil, slope, and aerial data (April flight), and secondly with crop, slope, and aerial data (July flight).
Fuzzy c-means was used in order to identify and further delineate in-field areas that possibly hold specific soil and crop characteristics and for which certain management strategies can be applied. Prior to the application of the algorithm, it was decided, in all cases (soil and crop characteristics), to preserve and output two zones per field for feasibility reasons in further managing them. Site-specific fertilizer management in order to be realistic and applicable should be adapted to farming conditions and to producers in terms of cost and labor effectiveness. However, except for practicality reasons, the final decision about the number of the management zones was also subject and yet controlled by the resulting between-zone variability or homogeneity in the analysis parameters. Sustainability 2021, 13, x FOR PEER REVIEW 10 of 25 Fuzzy c-means was used in order to identify and further delineate in-field areas that possibly hold specific soil and crop characteristics and for which certain management strategies can be applied. Prior to the application of the algorithm, it was decided, in all cases (soil and crop characteristics), to preserve and output two zones per field for feasibility reasons in further managing them. Site-specific fertilizer management in order to be realistic and applicable should be adapted to farming conditions and to producers in terms of cost and labor effectiveness. However, except for practicality reasons, the final decision about the number of the management zones was also subject and yet controlled by the resulting between-zone variability or homogeneity in the analysis parameters.

Spatial Variability Analysis
Towards investigating and capturing the spatial variability in soil and plant parameters for both fields, semivariograms were calculated as part of the ordinary kriging interpolation method. The best-fitted models were identified, describing the spatial patterns of different properties. In Table 5, the results of semivariogram analysis are presented. According to Table 5, the best-fitted models in the studied properties differed from one field to another in five of eleven properties (EC, Organic Matter, Total Nitrogen in soil, Potassium and Copper). In the other six parameters (Clay, pH, CaCO3, Total Nitrogen in plant, Iron and Zinc) the models best fitted between the two fields were the

Spatial Variability Analysis
Towards investigating and capturing the spatial variability in soil and plant parameters for both fields, semivariograms were calculated as part of the ordinary kriging interpolation method. The best-fitted models were identified, describing the spatial patterns of different properties. In Table 5, the results of semivariogram analysis are presented. According to Table 5, the best-fitted models in the studied properties differed from one field to another in five of eleven properties (EC, Organic Matter, Total Nitrogen in soil, Potassium and Copper). In the other six parameters (Clay, pH, CaCO 3 , Total Nitrogen in plant, Iron and Zinc) the models best fitted between the two fields were the same. The models selected, varied between spherical, pentaspherical, Gaussian, circular and exponential, reflecting the different spatial variability caused by the type of the property and also related to field conditions. A generic rule of thumb on which model is best linked to a specific soil or plant property, cannot be drawn from present analysis.
Based on the criterion that the ratio of nugget to sill is a measure of spatial dependency of the variables, according to Table 5, a large part of the variance in all parameters is introduced spatially. More specifically, as stated in [54], a ratio below 0.25 indicates strong spatial dependence in the variable. Following, the variance is medium if the ratio is between 0.25 and 0.75 and for ratio's values above 0.75 the spatial dependence is weak. Clay content, EC, CaCO 3 , Total nitrogen in soil and in plants, Potassium, Iron in Field 1 and CaCO 3 , Total nitrogen in soil, Potassium and Iron in Field 2, showed a strong spatial dependence. Medium spatial variance was revealed among pH, Organic Matter, Zinc in Field 1 and Clay content, EC, pH, Organic Matter, Total nitrogen in plants and Zinc in Field 2. Finally, weak spatial dependence was attributed only to Copper concentration in both fields. The spatially interpolated maps of Clay content, EC, CaCO 3 , Organic Matter, Total nitrogen in soil and in plants, Potassium and Iron in both fields exhibiting moderate to high spatial dependence are presented in Figure 8.   As far as soil properties are concerned, the northwestern part of Field 1 presents high values in clay percentage, carbonates, organic matter and total nitrogen, which can be attributed to its finer texture, while the southeastern area appears to retain lower values in these attributes. Electrical conductivity at the northwestern part of Field 1 presents some of the lowest estimated values, however low to medium values of EC are also widespread across the area of the field. On the other hand, in Field 2, electrical conductivity variability clearly follows the spatial structure of soil texture. Carbonates in Field 2, diminish from northeast to south, while organic matter and total nitrogen set a low value area, approximately at the middle of the field. Regarding plant properties, in both fields, total nitrogen and potassium concentration present only few scattered low spot areas, while the range between minima and maxima is not quite notable. Iron concentration in leaf samples on both fields however, reveals distinct spatial patterns. In Field 1, low concentration areas appear at northwestern, center and southeastern parts of the field and only two areas of high values appear mainly at the northern and western segments of the field. In case of Field 2, it is almost divided into two areas, based on iron concentration, one at northeast holding low values, and another at south to southeast having higher iron values.
Geospatial analysis of each studied parameter revealed valuable information on the way each variable fluctuates across space, and can serve in guiding soil and leaf sampling schemes. Moreover, the combined effect of all parameters along with near-infrared As far as soil properties are concerned, the northwestern part of Field 1 presents high values in clay percentage, carbonates, organic matter and total nitrogen, which can be attributed to its finer texture, while the southeastern area appears to retain lower values in these attributes. Electrical conductivity at the northwestern part of Field 1 presents some of the lowest estimated values, however low to medium values of EC are also widespread across the area of the field. On the other hand, in Field 2, electrical conductivity variability clearly follows the spatial structure of soil texture. Carbonates in Field 2, diminish from northeast to south, while organic matter and total nitrogen set a low value area, approximately at the middle of the field. Regarding plant properties, in both fields, total nitrogen and potassium concentration present only few scattered low spot areas, while the range between minima and maxima is not quite notable. Iron concentration in leaf samples on both fields however, reveals distinct spatial patterns. In Field 1, low concentration areas appear at northwestern, center and southeastern parts of the field and only two areas of high values appear mainly at the northern and western segments of the field. In case of Field 2, it is almost divided into two areas, based on iron concentration, one at northeast holding low values, and another at south to southeast having higher iron values.
Geospatial analysis of each studied parameter revealed valuable information on the way each variable fluctuates across space, and can serve in guiding soil and leaf sampling schemes. Moreover, the combined effect of all parameters along with near-infrared spectral data and slope variable were investigated and realized through PCA and fuzzy clustering methodology as presented in Figure 7.

Soil Characteristics
The degree of correlation between soil parameters, slope, and NIR as participated in the analysis was defined through Pearson coefficient (r) calculation. The correlation matrix is presented in Table 6. The majority of the parameters were correlated with each other in both fields. In Field 1, NIR was significantly and negatively correlated with clay content, organic matter, carbonates, and total soil nitrogen, while it was positively correlated with electrical conductivity. In Field 2, NIR shown significant negative correlation only with total soil nitrogen. On the other hand, increases in organic matter, carbonates, and soil acidity in Field 2 seem to increase NIR digital values. This diverse effect between the two fields may be attributed to noise in NIR values coming from existing weed patches and stones in the surface of Field 2. Among soil variables, electrical conductivity was negatively correlated with pH, organic matter, carbonates, total soil nitrogen in both fields, and with clay content in Field 1, which was also noted in spatial variability analysis (Figure 8). In Field 2, clay content was positively correlated with electrical conductivity and organic matter with carbonates and total soil nitrogen. Slope parameter was also correlated with all soil variables except for organic matter. Correlations among soil parameters, slope, and NIR indicated that PCA could be performed to aggregate those variables that set the main sources of variability across the fields.
In the case of Field 1, the principal components selected, were the first two according to the Kaiser criterion. Their eigenvalue was higher than 1 and the cumulative variance explained by those two components was higher than 72% of the total variance in the initial parameters (Table 7). More specifically, the first component accounts for 57.6% of total variance while the second one adds further 15.3% in the variance explained by the components extracted.  Table 8 presents each principal component's loadings. According to Table 8, the first component aggregates near-infrared (NIR) parameter with all soil parameters except for soil acidity (pH), which is highly correlated with the second component axis. This means, that NIR value monitoring and recording may infer crucial results regarding spatial variability in soil characteristics. It is also noticed that the values of NIR along with the conductivity of soil samples are grouped in the negative part of first component's axis, having the opposite influence on the component than that of the other soil parameters. This practically means that under the conditions of the first pilot field the variance of NIR values across the field is followed by the variance in electrical conductivity and inversely relates to the variability of the other five soil parameters, namely clay content, pH, CaCO3, OM, and total nitrogen.  Table 8 findings, in essence, constitute two linear combinations of the initial parameters with the resulting principal components, where the loadings represent the equations' coefficients. Subsequently, "solving" each equation for every single pixel, holding all seven values participating in the analysis, a pair of "scores" per pixel is calculated. Those pairs were introduced in Fuzzy c-means algorithm thus ascribing each pixel to a certain category according to its attributed, by the algorithm, membership value. The resulting layout, as produced in ArcGIS, is presented in Figure 9. In Figure 9 and according to the analysis applied in soil and reflectance parameters, Field 1 is divided into two zones (A and B). Zone A preserves higher percentages in organic matter, clay content and carbonates comparing to Zone B. Electrical conductivity is lower in Zone A than in B while total nitrogen content is slightly higher and soil acidity's difference is negligible. Mean values of all soil properties except for total nitrogen and pH were statistically different (p < 0.01) between the two zones (Table 9). Comparing Zones A and B, based on Table 9 results, Zone A presents better growth conditions for cotton plants as it can promote nutrient and moisture availability.   Figure 9. Zones recognized by Fuzzy c-means algorithm in Field 1 concerning soil characteristics. Table 9. Mean values of soil parameters per zone in Field 1. Regarding Field 2, PCA extracted in a total of eight components, of which the first three were selected representing 87.3% of the variance in the initial parameters, meaning that, counting only three variables, i.e., the three components, the variability of the eight initial factors can be fairly explained. The first component represents 41.4% of the variance in the participating variables, the second one adds 29.3%, and the third one contributes another 16.5% (Table 10). Consequently, the loadings of the two principal components for Field 2 are presented in Table 11. In the case of Field 2, the parameter of NIR is classified into the second component, forming a group with organic matter and carbonates. As previously mentioned, these two soil chemical parameters appear to be aggregated in the analysis of Field 1, along with the clay variable, which, in Field 2, is clearly grouped with electrical conductivity under the first component. The latter in soil substrate can be explained by the fact that heavy soils, rich in clay content, appear to have bad drainage conditions, which can cause the accumulation of salts and, thus, raise soil electrical conductivity. In contrary to Field 1, under Field 2 conditions, NIR parameter seems to drive positively the second component similarly to organic matter and carbonate content. The introduction of slope variable in the analysis may have caused this differentiation, which does not retain clear participation in any of the three components. Total nitrogen has excluded itself from the first two components, designating the third one, denoting that, under the conditions of this analysis, does not associate with the other soil variables. The respective mapping layout of the fuzzy clustering procedure is presented in Figure 10. The respective mapping layout of the fuzzy clustering procedure is presented in Figure 10. In Table 12, the mean values of the seven participating parameters per extracted zone in Field 2 are presented. Statistically (p < 0.01), the mean values of clay content, electrical conductivity, organic matter, carbonates, and slope were significantly different. Zone B, according to Table 12, had significantly higher mean value in carbonates and slightly higher in organic matter than Zone A. The other parameters' mean values did not appear In Table 12, the mean values of the seven participating parameters per extracted zone in Field 2 are presented. Statistically (p < 0.01), the mean values of clay content, electrical conductivity, organic matter, carbonates, and slope were significantly different. Zone B, according to Table 12, had significantly higher mean value in carbonates and slightly higher in organic matter than Zone A. The other parameters' mean values did not appear to differ among the two parts of the field; however, the slope in Zone B was greater than in Zone A. Soil conditions in Field 2, either in Zone A or in B, set no limitation for cotton growth. However, the resulting delineation of the zones can serve as a guide for future soil and fertilizing management practices. Further, as outputted by the statistical means between zones in both fields, the two-zone approach, may offer an acceptable trade-off between practicality and variability highlighting.  Table 13 presents the degree of correlation among plant parameters, slope, and NIR digital values in Fields 1 and 2. NIR demonstrated weak correlation with all parameters; however, plant variables among each other exhibited significant correlations. Notable is the negative relation of total plant nitrogen with potassium and copper and the positive one with iron and slope in Field 2. Further, copper was positively correlated with potassium in Field 1 and with zinc in Field 2. Potassium in cotton leaves was also negatively correlated with slope percentage in Field 2. Despite the statistically insignificant correlation of NIR with crop variables, PCA should be applied in order to reveal groupings and variances among data. The application of PCA in regards to parameters derived from plant tissue analysis in conjunction with NIR values taken by air, resulted in data presented in Table 14. Six components were calculated, among which, the first three cumulatively explain almost 69% of the variability in the initial data. This is a rather weak percentage judging by the fact that the parameters involved were only six and the selected components were three, thus reducing the dimensions of the analysis only by 50%. However, the first two components can explain 51.5% of the variance in the parameters engaged. Furthermore, a few groupings were revealed, as shown in Table 15. In fact, potassium and copper concentrations in cotton plants were associated with the first principal component, total nitrogen and zinc concentration had high correlation with the second one, and NIR values, followed by iron concentration were grouped under the third component. It is clear that, under the specific dataset, the NIR parameter cannot be used for drawing safe conclusions about the spatial distribution of the nutrient composition of the plants throughout Field 1. The mapping of the fuzzy classification led to the layout presented in Figure 11. Contrary to the zones extracted by the soil characteristic analysis (Figure 9), here the zones are not aggregated in specific parts of the field, but the sub-areas are scattered across the field. It is also noticeable that the algorithm even classified cotton rows into the two zones. For this reason, the proposition of site-specific management strategies regarding crop nutrition is rather difficult to be applied.   Figure 11. Zones recognized by Fuzzy c-means algorithm in Field 1 concerning crop characteristics. Figure 11. Zones recognized by Fuzzy c-means algorithm in Field 1 concerning crop characteristics.

Crop Characteristics
Further, it is argued that, according to the results presented in Table 16, clustering did not manage to discriminate areas of clear variability between the crop parameters measured, except for iron concentration, in which statistically significant differences (p < 0.01) in mean values per zone were recorded. Zone B appears to retain higher mean values of iron comparing to Zone A, however, being at absolute values that set no risk of future deficiencies in the plants. In the case of Field 2, with the inclusion of slope percentage, the analysis extracted seven components from which the first three accounted for 72.6% of the variance in the data (Table 17). The NIR parameter seems to have low contribution to all of the three components; however, according to its loadings, tends to group under the first component. This states its slight correlation with the other parameters measured. The variables of slope and potassium concentration are also grouped under the first component having high absolute loadings. The loadings of zinc and copper clearly differentiate from the other parameters, characterizing the second component, while the third component is dominated by an iron parameter and, at a lesser degree, by total nitrogen (Table 18). According to the results of PCA in crop characteristics, NIR did not group in sufficient rate with the other parameters in Field 1 and in Field 2. A possible explanation is the fact that all nutrients measured in cotton plants in both fields were at sufficient levels and their value variations along the area of the fields were not significant. Due to this uniformity, the variability of aerial NIR values was also low, which was reflected in the output of PCA. The mapping layout of the Fuzzy c-means output regarding Field 2 is shown in Figure 12. Visually, there is a clear delineation of two zones, one at the left-hand side and the other at the right of the field. The bottom right corner of the field being classified in Zone B is possibly attributed to the shade effect that a nearby tree had on the background reflectance, and it is considered an image artifact. Zone B is possibly attributed to the shade effect that a nearby tree had on the background reflectance, and it is considered an image artifact.  (Table 19) did not statistically differentiate among the two zones, with the exception of iron concentration (significance level p < 0.01), as similarly observed in Field 1. In both cases, despite the minimum variability observed between the two zones extracted per field, zoning-as achieved-may contribute and guide future plant sampling schemes.   (Table 19) did not statistically differentiate among the two zones, with the exception of iron concentration (significance level p < 0.01), as similarly observed in Field 1. In both cases, despite the minimum variability observed between the two zones extracted per field, zoning-as achieved-may contribute and guide future plant sampling schemes.

Conclusions
In the context of this study, UAS data were analyzed with soil and crop parameters in two cotton fields during a growing period. All data were analyzed using geostatistical and geospatial methodologies coupled with PCA and Fuzzy c-means clustering under GIS environment. A primary conclusion, of practical importance is that the use of UAS in crop production offers a quick and reliable way to monitor soil and plant capital. As many UAS manufacturers throw their products into the market, raising competition, and diminishing costs, many farmers are able to engage UAS or similar services to their production procedure. Another finding is that fieldwork, in the form of sampling and further analyzing soil and tissue samples, are of crucial importance for evaluating and confirming remotely sensed data. This work demonstrated that reflectance data were correlated with organic matter, carbonate, and clay content, while those data can be directly used for in-field zone delineation. This time efficient form of field monitoring can lead to targeted sampling schemes, avoiding unnecessary work and costs. On the contrary, crop nutrient characteristics were insignificantly combined with aerial data, which may have been the result of uniformity of the measured parameters across the two fields. The resulting zones, as mapped in this case, offer little in fertilizing management guidance. Further work should evaluate present results by incorporating more fields, crops, and growing seasons. Research efforts should also focus on translating the outcomes of soil and crop monitoring through expert decision-making tools and on efficiently applying management plans through variable rate technologies. The integration of new technologies and established methodologies in primary production, such as those demonstrated, provide notable means for applying site-specific crop management in broader adoption levels and constitute a crucial motive for visualizing, designing, implementing, and assessing environmental strategic plans towards a circular economy.