QVigourMap: A GIS Open Source Application for the Creation of Canopy Vigour Maps

In a precision agriculture context, the amount of geospatial data available can be difficult to interpret in order to understand the crop variability within a given terrain parcel, raising the need for specific tools for data processing and analysis. This is the case for data acquired from Unmanned Aerial Vehicles (UAV), in which the high spatial resolution along with data from several spectral wavelengths makes data interpretation a complex process regarding vegetation monitoring. Vegetation Indices (VIs) are usually computed, helping in the vegetation monitoring process. However, a crop plot is generally composed of several non-crop elements, which can bias the data analysis and interpretation. By discarding non-crop data, it is possible to compute the vigour distribution for a specific crop within the area under analysis. This article presents QVigourMaps, a new open source application developed to generate useful outputs for precision agriculture purposes. The application was developed in the form of a QGIS plugin, allowing the creation of vigour maps, vegetation distribution maps and prescription maps based on the combination of different VIs and height information. Multi-temporal data from a vineyard plot and a maize field were used as case studies in order to demonstrate the potential and effectiveness of the QVigourMaps tool. The presented application can contribute to making the right management decisions by providing indicators of crop variability, and the outcomes can be used in the field to apply site-specific treatments according to the levels of vigour.


Introduction
According to the definition adopted by the International Society of Precision Agriculture, Precision Agriculture (PA) can be defined as the "management strategy that gathers, processes and analyses temporal, spatial and individual data and combines it with other information to support management decisions according to estimated variability for improved resource use efficiency, productivity, quality, profitability and sustainability of agricultural production". [1] Sassu et al. [2] argue that PA is related with the use of technology to manage the spatial and temporal variability associated with agricultural production, so that it can improve the crop performance, bringing economic benefits and environmental quality with the correct use of pollutants. Regardless of the sources used to establish PA, they all include several innovative technologies, such as the use of data acquired by Unmanned Aerial Vehicles (UAVs), satellite imagery data, the use of Geographical Information Systems (GIS) technologies, nutrient management field mapping, Internet of Things (IoT) sensors, and data processing techniques such as Machine Learning and Deep Learning, among others, allowing experts to use specific described. The experimental results of the case studies and performance evaluation are presented in Section 3. The discussion and conclusions are addressed in Section 4.

Methodology
The developed application relies on a methodology with different steps (Figure 1). Each step is responsible for the generation of a different outcome, or for the performance of different data operations. Raster products resulting from the photogrammetric processing of UAV multispectral data can be used as the input. These are the cases of radiometrically corrected spectral reflectance or VIs-of which the present crop status indicators are used-and Crop Surface Models (CSMs), from which geometrical properties can be extracted. Within a crop field, there are other elements beyond crops, such as bare soil, shadows, and invasive and/or undergrowth vegetation [28]. The inclusion of information from those elements can bias the results, raising the need to mask them out. In order to discard information from non-crop elements, a vegetation segmentation step needs to be performed, preserving the majority of the information related to the crop to be analysed. Among the available segmentation techniques, image thresholding is one of the simpler methods [29]. In this way, a given threshold value can be used to transform a single-band image into foreground and background elements, in which the values from the foreground are kept and the background elements are discarded from the further steps [28]. This process is represented in Equation (1), where the value of a given image A at the pixel in the row i and column j is verified with regard to whether it is greater or equal than the threshold value T, creating a new image M with the pixel values from A when those values are greater or equal than T, or with no data values (NaN) for pixels below T. , , , This process is applied to all of the pixels of a single-band raster, and can be implemented in multiple UAV-based outcomes. The threshold value T can be inferred from the histogram analysis of the study area or by a predefined value, because in most cases the values from VIs tend to be higher in crops than in other vegetation or bare soil. If a geometrical outcome is being used, as in the case of CSM, T can represent a minimum height range, discarding all of the pixels that do not satisfy this condition. The combination of multiple outcomes can also refine the crop segmentation step, a raster with only the values of a VI within a certain threshold and a minimum height. If a height threshold is used, other vegetation (undergrowth) can be identified with pixels that were not within the minimum height but satisfied the threshold value of the VI.
After the segmentation stage, the image M is subjected to a filtering process in which the values are interpolated and then smoothed, filling the gaps left as no data by the segmentation stage. Before the computation of the vigour map, vigour classes can be used to reclassify the values, which can be performed manually, according to a given interval, or automatically, according to the data distribution in the quartiles. Thus, a vigour map is produced with a given number of classes which best represent the study area.
Apart from variability mapping, vigour maps can also be used to compute other outcomes, such as prescription maps to be deployed on the field assisting in the application of treatments. This can be carried out by assigning each vigour class to a rate to be applied in the field.
Moreover, other important outputs can be obtained, as in the case of the vegetation distribution map. If height information is available, this map can be distributed in three classes: one representing the vegetation under analysis, based on the pixels that meet the threshold values of both VI and height; another representing other vegetation in the parcel, mostly undergrowth, which satisfies the VI threshold but does not have enough height; and one class with areas without vegetation. If height information is not used, two classes can be obtained, one representing vegetated areas and another with no vegetation parts.

QVigourMaps Application
In the last decade, FOSS has emerged in the GIS context. Several pieces of GIS software, such as the Geospatial Data Abstraction Library (GDAL/OGR) [30], Geographic Resources Analysis Support System (GRASS) [31], and System for Automated Geoscientific Analyses (SAGA) [32], have been updated with new functionalities over the years. QGIS software is one of the most used GIS, being user-friendly, easy to use, intuitive and it is complemented with algorithms from other software such as GRASS or SAGA. This complement is integrated from the Processing Toolbox, which is composed of a list of all of the available algorithms grouped in different blocks called Providers, and custom models and scripts that can be added in order to extend the set of tools [24]. Considering all of the advantages of using the QGIS environment, a GIS open-source application named QVigourMaps was created under QGIS software.
A set of functional requirements were defined for QVigourMap: (1) it should be able to produce multiple outcomes (vigour maps, vegetation distribution maps, prescription maps based on the vigour levels); (2) it must contain a GUI that enables it to select the inputs to be used and allow a certain level of parameterization for the refinement of the obtained results; and (3) it must be able to produce results from any provided VI, and the use of height information can be optional. As non-functional requirements, it was defined that the developed software should be open-source, time efficient, have low computational power requirements, and be able to be performed on an intermediate level laptop.
The QVigourMaps application/plugin was developed using the Python programming language [33]. In order to start the plugin development, the Plugin Builder was used. The code was developed using the PyCharm editor [34]. The main GUI was created as a dialog button composed of five sections: (i) Vegetation Index; (ii) Crop Surface Model; (iii) Study Zone, with two approaches, i.e, using a shapefile or using the Preview Study Zone window; (iv) Resampling Parameters; (v) Output section, including two different outputs, i.e., the final vigour map and the vegetation distribution map ( Figure 2). The classification map is generated by automatically reclassifying the vigour map in three different classes: 0-Soil, 1-Other vegetation, and 2-Crop.
The selection of QGIS software to integrate QVigourMaps was due to the fact that this GIS is open-source and has a broad community using it, along with the possibility to use its GIS environment to perform other tasks or to use the other tools and functions available.
The Vegetation Index section allows us to define the input directories of an alreadycalculated VI map, or to directly select a layer already opened in QGIS. The Crop Surface Model input section was implemented with the possibility to calculate the CSM providing a Digital Surface Model (DSM) and a Digital Terrain Model (DTM) as inputs or directly giving an already-computed CSM. The study zone can be defined through two options: by a shapefile or by using the Preview study zone window. This window displays the input maps inserted into the input sections. The user can then define the study zone using the cursor by drawing a polygon. An example of this procedure is displayed in Figure 3. The red polygon shown in Figure 3 represents an example of a study area. It was defined under the CSM map, but it can be defined using another input map. The application internally converts the polygon to a shapefile and clips all of the input data to that area. All of the intermediate results are saved in a temporary folder that is removed when the user closes the application.
The Resampling parameters section allows us to define the parameter values in all of the procedure, and it is composed of the Horizontal Filter Radius and the Vertical Filter Radius value; the minimum value of VI (Minimum value (Vegetation Index)) and the minimum height from CSM (Minimum height (CSM)). By default, the values assigned are 2 for the horizontal and vertical filter radius, and 0.5 for both minimum values, but can be modified by the user. It should be noted that the values used for the resampling filter radius must be selected according to the planting parameters of the crops, and according to the spatial resolution of the data. If the filter radius is too high, it will cause an excessive data smoothing that would prevent the effective analysis of the crop vigour variability.
On the other hand, if a low filter radius is used, the final vigour map would look similar to the raster used as the input of this step (interpolated vegetation mask with NDVI values), with no or minimal smoothing applied.
In the Output section, two raster outputs are defined: the vigour map, which is saved in a specified directory, and the output vegetation distribution map with three main classes, i.e., soil, other vegetation, and crops. After that, the user can, optionally, reclassify the results using the Reclassify map button. It opens a new window composed of: (1) a table with four columns; (2) a button to add more rows to the table; (3) a button to remove rows from the table; (4) a button to recalculate the quantile intervals if the number of rows is modified; (5) a textbox to define the output path of the reclassified vigour map and another for the prescription map ( Figure 4). The table cells are editable, and are composed of four columns for the values inserted in each row, including: (1) the minimum interval value; (2) the maximum interval value; (3) the vigour class value; and (4) the prescription value.
The minimum and maximum values are automatically defined in three quantile intervals. This procedure was implemented in the code using the numpy library and its percentile function. These values are estimated based on raster values and using GDAL functions. If the user modifies the number of rows (adding or removing rows), the button Recalculate quantile button must be used to estimate new quantile intervals. The vigour class values and the prescription values must be inserted manually by the user. In order to ascertain the processing time, a progress message bar was incorporated in the application. The progress bar was implemented using the progressMessageBar widget and the assignment of time value through time Python module.
One of the strong advantages of QGIS software is the integration of GIS algorithms from external GIS software, such as GRASS and SAGA. Several algorithms are available in the QGIS Processing Toolbox. For instance, grass7:r.map.calc.simple (from GRASS) and saga:cliprasterwithpolygon (from SAGA) allow us to limit the raster to a study zone (mask). Furthermore, grass7:r.resamp.filter (from GRASS) and saga:resamplingfilter (from SAGA) allow us to resample a raster. There are also native algorithms from QGIS software and GDAL algorithms, that have related functionalities, such as gdal:cliprasterbymasklayer and the clip algorithms from GRASS and SAGA. In this work, similar algorithms were manually tested in order to implement the algorithm with the best performance regarding the execution time and the parameters involved. Taking these aspects into account, in the application, different algorithms from GRASS, GDAL and SAGA (from QGIS Processing Toolbox) were used, such as grass7:r.map.calc.simple (to perform the clipping with the study area), saga:closegaps (to eliminate the nodata values) and grass7:r.resamp.filter (to resample a raster map using an analytic kernel).
The application is free, open-source and available on GitHub (https://github.com/liaduarte/QVigourMaps, in Supplementary Material), contributing also to the geoinformatics advance in PA solutions.

Case Study
For the purpose of testing QVigourMaps in different contexts, two study areas were selected: a maize field and a vineyard. These areas are located within the campus of the University of Trás-os-Montes e Alto Douro (Vila Real, Portugal), and are highlighted in Figure 5. The vineyard plot (approximately 0.6 ha) is composed of red grapevine varieties typical to the Douro Demarcated Region. The maize field (0.5 ha) is used for silage production. The selection of these areas is based on the different levels of vigour and missing plants, thereby enabling the developed application to be tested in a multitude of scenarios. Figure 5. Location of the study areas used for the testing of the developed application (a) and its overview (b). The orthophoto mosaic used in (a) is not representative of the data status used in this study.

Data Acquisition
The multispectral data used in this study were acquired using the multi-rotor DJI Phantom 4 (DJI, Shenzhen, China) with a Parrot Sequoia (Parrot SA, Paris, France) mounted in the UAV. This camera was composed of a 16 MP RGB and four 1.2 MP multispectral imaging sensors in order to acquire single-band images in the green, red, red-edge (RE), and near-infrared (NIR) parts of the electromagnetic spectrum. For the radiometric calibration, irradiance data is acquired during the flight, and reflectance data is acquired using a calibration target. For the purpose of this study, RGB imagery was not used. The data was acquired in seven flight campaigns: four in 2018 (maize) and three in 2019 (vineyard). Table 1 provides information regarding the specific date of each campaign and the flight parameters for each crop.

Data Processing
The multispectral data acquired in each flight were processed using Pix4DMapper Pro (Pix4D SA, Lausanne, Switzerland). A high-density point cloud was generated, and from its interpolation, a DSM and then a DTM were computed. After the performance of the radiometric calibration, reflectance maps of each band are generated, and the NDVI [11] is computed using the reflectance from the NIR and red bands, as in Equation (2).

NDVI
NIR Red NIR Red A CSM was computed from the subtraction of the altitude values of the DTM to the DSM, as in Equation (3).

QVigourMap Use Cases
In order to exemplify and evaluate the usage of the QVigourMap, the two case studies described in Section 2.3 were used. In the vineyard plot, the vegetation dynamics were monitored and the generated vigour maps were compared with the ones computed without applying the vegetation segmentation procedure. The case study of the maize field was used to monitor the crop status and vigour changes along its growth.

Vineyard
For the vineyard case study, none of the pixels with a value bellow 0.4 were considered, and a value of 0.3 m was used for the height threshold in the CSM. In this crop, the value 2 was used for the horizontal and vertical radius. These parameters were applied in all of the flight campaigns. Figure 6 shows the inputs used in QVigourMap using the data from the first flight camping. As can be seen in Figure 6c, the vegetation segmentation step enabled us to divide the grapevine vegetation from the other features present in the vineyard plot. Because there is a vegetation segmentation step in the QVigourMap application, it is possible to compute a vegetation distribution map that can be used to understand the vegetation growth/decline dynamics, and to identify areas with missing/little vegetation. Figure 7 presents the vegetation distribution along the three evaluated flight campaigns. Between the first flight campaign to the other two there are two main differences: the grapevine vegetation presents a small decline while other vegetation increased.
The vigour maps generated for each flight campaign are presented in Figure 8, using three levels according to quantile intervals (high, medium, and low). For comparison purposes, two vigour maps were generated for each flight campaign, one by directly resampling of all of the pixels from NDVI (top row of Figure 8), thereby including information from other vegetation and bare soil, and another using the approach applied by QVigourMap, considering only the grapevine vegetation (bottom row of Figure 8) higher than a minimum NDVI value and with a certain height. The resampling was achieved by using a filter with a box kernel, and with a vertical and horizontal radius of three.  If we observe the vigour maps generated with both approaches, there are common vigour areas. Low vigour is most noticeable in the bottom part of the analysed vineyard plot, medium vigour is noticeable towards the centre, and higher vigour is verified in both the northwester and eastern parts of the vineyard. By removing information related to soil and non-grapevine vegetation (filtered approach), we enabled the mapping of the real grapevine vigour levels, and the low and high vigour areas in the unfiltered approach presented different vigour levels in the filtered approach. Analysing the temporal vineyard variability, the expansion of the area of high vigour in the northwestern part of vineyard was verified from the first to the remaining flight campaigns. The high vigour area in the eastern part of the vineyard decreases along the flight campaigns. These effects were verified in both approaches, with more incidence on the filtered approach at the last two flight campaigns. The low vigour area in the southeastern part remains nearly stable in the unfiltered approach, while the filtered approach presents areas with medium and high vigour.

Silage Maize
Regarding the monitoring of the silage maize growth, only the NDVI was used as the input, with a different threshold value considered for each flight survey; this value was increased by 0.1, stating at 0.3 in the first and ranging up to 0.6 in the last campaign, before sowing. These values were selected according to the NDVI value of the vegetation (top row of Figure 9). A box filter with a vertical and horizontal radius of two was used for the data resampling. As expected, the overall NDVI value increased along the growth. A higher heterogeneity was observed in the first two flight campaigns, while this is less noticeable in the third and almost not perceptible in the final survey. However, despite the overall increase on the NDVI values and the apparent heterogeneity decrease, the vigour zones remain reasonably stable throughout the flight campaigns (bottom row in Figure 9). Generally, the eastern part of the plantation presents a lower vigour than the opposite site, with a medium vigour transition among them. Small variations are observed from the first to the second flight campaign, whilst the comparison of these campaigns with the last two verifies a settlement of the higher vigour area. From the third to the final flight campaign the low vigour areas decreased, transitioning to medium and even high vigour in some parts.

Processing Time Performance
Two of the goals of this QGIS plugin are to work with few parameterizations and to require low computer processing power in order to easily operate in different computer environments. Therefore, the performance of QVigourMap was evaluated by analysing the time spent in each step of the processing pipeline. For this purpose, the results from one campaign for vineyard (flight 1 DOY 206) and for silage maize (flight 3 DAS 52) were used. The procedures were performed using QGIS 3.18 Zurich installed on a Windows 10 operative system running in a computer with an Intel ® Core™ i7-7500U CPU @ 2.70 GHz 2.90 GHz CPU, 8 GB RAM (2133 MHz), and with a 256 GB SSD (HFS256G39TND-N210A). The time spent in a manual procedure, executing each individual operation until the same final result is reached, was also measured for further comparison with QVigourMaps. The results of the performance tests are presented in Table 2.
In the manual procedure, the following steps are required to be performed: (i) Clip the raster inputs with a polygon; this step is performed before the main procedures, because if we have larger input files, the restriction to the study zone will minimize the time procedure. (ii) Both the vegetation mask and height mask are similar procedures, which use the r.mapcalc.simple algorithm, in which each of them allows us to create a mask of the VI and CSM in a specific interval of values (these values are chosen by the user in the GIS application), as in (1), in which each pixel (i, j) assumes the value of "1" or "NaN". (iii) The intersection of vegetation and height masks refers to the combination of the CSM binary mask (H)-if used-and the VI mask (V) to produce a new raster (N) containing vegetation index values (VI) with the pixels with a value of "1" in both masks, as in (4), this operation is performed using the r.mapcalc.simple algorithm. (iv) The vegetation interpolation is performed by using the close gaps algorithm from SAGA, which was used to eliminate the no data (Nan) values in the result obtained. (v) Data resampling by using the resampling filter algorithm from GRASS allows us to resample a raster map using an analytic kernel.
, , , The total time spent using the manual procedures, i.e., performing a step-by-step approach using algorithms from the Processing Toolbox in QGIS was 3′55.16″ (235.16″) in the vineyard dataset and 2′48.88″ (168.88″) in the maize dataset. The masking operations and closing the gaps took most of the time in both datasets. These performance measurements include the estimated time for the user to insert the inputs, outputs and parameters in each tool. The total time spent to generate a vigour map in the vineyard dataset using QVigourMaps was 1′45.86″ (105.86″), this being the difference between the two procedures of 2′09.30″ (129.3″). In the maize dataset, the QVigourMaps application created the vigour map in 1′40.25″ (100.25″), for a difference of 1′08.30″ (68.30″). It can be noted that the total time presented in the vineyard and in maize crops includes not only the vigour map creation but also the vegetation distribution in the case in which the CSM map was introduced (vineyard). Table 2. Comparison of the performance of QVigourMaps through the processing times with the manual procedures.

Crop
Step

Discussion and Conclusions
In this work, an open-source application, QVigourMap, developed under QGIS software was presented. The application was implemented using algorithms from the Processing Toolbox (GRASS and SAGA). This is a great advantage of using QGIS software. Its other advantages are that: (i) it is free to use; (ii) it is intuitive, as is presents a GUI which is easy to understand; (iii) it minimizes the time to process all of the steps (if they are performed step by step); (iv) it has a tutorial to support the user; (v) it can be updated at any time, with other functionalities, and by any other user/developer; and (vi) the code is hosted in GitHub.
The QVigourMap application and its integration in QGIS as a plugin opens the possibility for an easier interpretation of VIs through vigour maps by filtering nonvegetation elements. The potential of this tool for multi-temporal vegetation monitoring was explored in the vineyard and maize case studies.
The results presented in Section 3 demonstrate the potential of the QVigourMap to generate vigour maps by masking out other vegetation and bare soil, as seen in Figure 6c. Regarding the vineyard case, the vegetation changes through the period analysed ( Figure  7) were related with the grapevine phenology stages, because the data covers the whole period from the beginning of veraison to harvest time [35]. Moreover, the occurrence of precipitation in the days prior to the data acquisition justified the growth verified in the other vegetation in the last two flights [29]. Vegetation distribution maps can be used for soil corrections and/or to plant/replace crops, as in the vineyard case.
Because non-grapevine information was discarded, the vigour maps produced with the approach used in QVigourMap, when compared to an unfiltered approach (Figure 8), presents a better variability discrimination within the vineyard. Other studies demonstrated that the vineyard vigour correlates with crop height [36] and canopy temperature [12]. Moreover, the vigour maps generated by considering the whole information can still provide an overall perspective of the vineyard. Some studies explored different vigour areas in multi-temporal datasets, with vigour areas are defined only by analysing the vineyard in a single flight [37]. As demonstrated in this study, vineyard variability can radically change during the same season; this trend was also verified in data from different years in different seasons [38].
With regard to the maize plantation (Figure 9), different NDVI values were considered for vegetation filtering; these values were selected in order to consider the majority of the crops in the field because it is well known that NDVI tends to present higher values in areas with higher healthy vegetation density [39]. In this case, CSMs were not used due to the nature of the performed flights and the resolution of the multispectral sensor, which resulted in a lack of spatial resolution to infer height variations from the plants at early stages. Thus, values from some weed patches could take up part of the vigour map. The vigour levels observed are nearly similar in the first two flight campaigns; this can be explained by considering the small temporal difference between both campaigns (seven days) and the early development stages. In the acquired data, after a period of greater growth-from the second to the third flight campaign (17 days)-the plants were already notorious. However, there are differences in the low vigour areas, in the right side of the plantation, which can be justified by a soil nutritional deficit. From the results obtained in the last two flights, there is a greater uniformity in the plantation; still, it is clear that on the right side of the plantation there are signs of stress, clearly influencing the plant vegetative status in that area. It should be noted that an irrigation system was implemented in this parcel, and that it was already in operation at the time of the last two flights [40]. Furthermore, nitrogen corrections were carried out, which help to justify the vigour changes in the low vigour area in the last flight campaign. Despite the changes verified in the last flight campaign, the overall vigour areas were maintained (i.e., the western part showed higher vigour than the eastern part of the plantation), which lead us to conclude that, in an early stage, the areas with higher production could already being noted. Thus, it can be concluded that, using multispectral data acquired by UAV, it is possible to check crop variability and to optimize the treatments or irrigation, as a measure of correction/mitigation of the culture under analysis.
The possibility of creating prescription maps in a vector format was not demonstrated in the presented case studies due to possibility to apply the treatments in the field in an automatic way. However, this functionality was implemented in QVigourMap by using the generated vigour map and resampling the vigour level to a certain rate for field application. This type of map can be used in UAV spraying applications, because UAVs can access other areas that are not reached by ground vehicles [41][42][43]. Moreover, the vigour maps presented in this study only used three vigour levels (high, medium and low) distributed in quantiles, and there is still a possibility to modify the number of vigour classes and to select the intervals for each class, as used in Pádua et al. [12]. As such, QVigourMap can be an excellent tool for decision support system agricultural management purposes.
Despite the fact that in the presented case studies only the NDVI was used, the application supports other VIs, which can be used depending on the desired output. Some Vis, which use reflectance from other spectral regions, can present a higher agreement for the estimation of certain biophysical parameters (e.g., leaf pigments or nutrient estimation) and the detection of crop variability [44][45][46]. Additionally, remote sensed data from spaceborne or airborne platforms can be used in the developed application for studies at larger scales.
Several conclusions can be drawn by analysing the results presented in Table 2. In the processing time evaluation, we can conclude that the GIS application is effective and not time-consuming. When compared with a manual approach, gains of 45% and 59% were reached in the vineyard and maize datasets, respectively. Besides this, a user with no experience with QGIS tools, especially the Processing Toolbox algorithms, would take a considerable amount of time to reach the same results. In this way, it can be affirmed that this application is more intuitive and easier to use, with a great potential of application for different crops. It can also be easily adapted to another methodology or use in other countries. The code is open and can be modified according to the user requirements. This is the basic idea of the open source.
In the future, more functionalities will be added to the QVigourMaps application, such as: (i) the incorporation of a dynamic temporal variation of the vigour maps obtained; (ii) improvements in the methodology; (iii) improvements in the application's usability, providing more customization options to the user; (iv) the implementation of quantitative vigour indicators (such as leaf area index); (v) the integration of new analysis tools. Moreover, a case demonstrating the potential of variable rate applications can be addressed to fully explore the potential of QVigourMaps from a multi-temporal perspective for full crop monitoring, in order to provide decision support insights and to assist in field operations.