Estimating Carrot Gross Primary Production Using UAV-Based Multispectral Imagery

Abstract

Gross primary productivity (GPP) is an essential parameter to assess the efficiency of terrestrial ecosystems on carbon transfer. Although GPP is regularly measured with eddy covariance (EC) systems, these are restricted to the tower footprint area, and remote sensing (RS) products have estimated GPP using multispectral vegetation indexes (VIs) from farms to whole ecosystems. Indeed, nowadays, unmanned aerial vehicle (UAV)-based RS technology is becoming more accessible. Accordingly, we propose the estimation of GPP using VIs at high spatial resolutions using UAVs and multi-spectral cameras. A small typical farm in Colombia was cultivated with carrot as our base crop. An EC system was installed to estimate GPP and was used as a reference. A total of 24 VIs from UAV-based RS products were selected and compared with the GPP of the EC system. A cross-validation process was performed, and seven VIs obtained a high R² score (0.76–0.78). The accumulated GPP, estimated with the best index (NIRv) was 520.3 g C m⁻², while the GPP-EC estimate was 580.4 g C m⁻² (10.3% error). This work showed that it is possible to estimate the GPP of carrot crops using UAV-based RS, VIs, and linear regression models, which can be used in further research on GPP using UAVs.

1. Introduction

Understanding the global carbon cycle in the context of climate change is currently a key point for researchers and society, and a critical parameter related to the carbon cycle is the terrestrial gross primary productivity (GPP). The accumulated GPP is used as a variable to calculate, over a period of time, the total amount of carbon dioxide fixed through photosynthesis [1,2]. The GPP is also used to assess the efficiency of terrestrial ecosystems on carbon transfer [3,4]. Changes in the GPP (diurnal and seasonal) depend on nutrient availability and climatic conditions (precipitation, humidity, light, and temperature), and the latter also influences its spatial distribution [5]. Estimations of GPP are decisive in monitoring carbon exchange in forest lands and agricultural ecosystems for the following reasons: information about GPP is essential for measuring and checking the variations in productivity; estimations have a spatial and temporal distribution; and the products can be used in optimizing agricultural management practices, especially over a global climate change context [6,7]. The GPP can be derived using the eddy covariance (EC) method as the difference between the net ecosystem exchange (NEE) and the ecosystem respiration (R_eco) during daylight [8]. Additionally, GPP can be obtained by using flux chambers; however, they are rarely used due to their small coverage and inaccurate results, associated with the collar insertion methodology [9].

The high spatial variability of the GPP affects its measurement from the ground at regional scales, making it very challenging [10]. For instance, the spatial coverage of the GPP estimated using the EC system is delimited to the tower footprint area [2]. Nevertheless, remote sensing (RS) has been used to estimate GPP from satellite platforms at different resolutions, contributing information to decisionmakers for improved management [2,11]. Spatially-continuous GPP estimates have been mainly obtained from three RS models [12]. First, there are the vegetation index statistical (VI) models (VI-GPP) that report a direct relationship between vegetation indexes (VIs) and GPP [10,13,14,15,16,17,18,19], as well as VIs combined with atmospheric variables, such as the photosynthetically active radiation (PAR) [1,18,20], land surface temperature [21], the land surface water index (LSWI), and the enhanced vegetation index (EVI) [22], among other models [12]. Second, there are light use efficiency (LUE-GPP) models, where the GPP is computed as the production of absorbed photosynthetically active radiation (APAR) by the plant canopy and the efficiency of converting APAR into carbon fixed trough photosynthesis, such as the MODIS-GPP model [23], the vegetation photosynthesis model (VPM) [10,16], the CI-LUE model [24], and several more LUE-based models [12]. Third, there are process-based models [5]. These models use complex model structures with climatic, biological, and soil variables to explain the main components of ecological processes related to photosynthesis [15,25]. The VI-based models are typical in agriculture because they reflect the actual status and growth stage of crops [26]. Estimation of GPP from VI-based models has been achieved through statistical relationships, using simple algorithms and few input variables [12]. Several researchers [2,5,10,13,14,15,16,17,18,19,20,23] have shown that VI from RS technology could be used to predict the GPP.

In this context, VIs are spectral imaging transformations (mathematical expressions) of two or more bands. These VIs are usually calculated as ratios, subtractions, or linear combinations of the light reflectance in the visible (VIS), near-infrared (NIR), and shortwave infrared (SWIR) spectral regions [27]. The normalized difference vegetation index (NDVI), enhanced vegetation index (EVI, EVI2), and near-infrared reflectance of vegetation (NIRv) are the most common indexes used to estimate GPP [28]. Other, less commonly used VIs correlating with GPP have been also proposed, including CIgreen [17,18], CIred-edge [20], red-edge NDVI (RENDVI) [20], the modified chlorophyll absorption ratio index (MCARI) [20], MERIS terrestrial coefficient index (MTCI) [20], soil adjusted vegetation index (SAVI) [1], water difference vegetation index (WDVI) [1], wide dynamic ratio vegetation index (WDRVI) [17], WDRVI green (WDRVIgreen) [18], land surface water index (LSWI) [10,22], and green–red vegetation index (GRVI) [29]. The spectral data used in the computation of VIs can be the result of different types of data acquisition platforms, such as spaceborne (satellite sensors), manned airborne vehicles, or unmanned aerial vehicles (UAVs) [26]. The resolution of satellites is relatively low, and the possibility of cloud cover limits the retrieval of primary productivity data at the field scale [30]. The case of manned airborne-based RS technology is limited by its temporal and spatial resolution and high cost [31]. In recent years, there have been important developments in UAV-based RS technology (or drones) in terms of cost, weight reduction, and operability (e.g., pre-programmed flights) [26,32,33]. As such, UAV-based RS can increase the scale and the resolution of primary productivity observations, overcoming the shortcomings of satellite and manned airborne platforms [30]. It is also important to have ground truth data, since different vegetation species will have different contributions to carbon fluxes, depending on the time of the year [9].

While the ability to use Vis, such as NDVI, EVI, and NIRv to track GPP using satellite-based platforms has been well established [28], the operational capability to accurately estimate GPP at the farm scale using UAV platforms has not been widely analyzed [30]. In particular, GPP estimates from remote sensing platforms and VIs have been conducted for maize [1,10,15,18], paddy rice [16], soybean [15,18], wheat [20], evergreen and deciduous forests [14,21], woody savanna [21], broad-leaved pine mixed forest, temperate steppe, alpine shrubland, alpine marsh, and alpine meadow steppe [22]. They found linear relationships between the daily GPP and VIs over different study areas. Less attention has been paid to deriving GPP from spectral data for vegetable crops using UAVs. The carrot crop (Daucus carota L.) plays a very important role in food security in Colombia and is one of the traditional consumption crops, with around 11,000 hectares planted in 2017, mainly by small producers, with an average yield equal to 30 t ha⁻¹ [34]. The main producing regions are Cundinamarca, Boyacá, Antioquia, and Nariño.

In this context, the present study aims to estimate the GPP of carrot crops from UAV-based multispectral VIs. The performance of this relationship has been evaluated based on field data collected during one growing season (2021) through the EC system measurements.

2. Materials and Methods

2.1. Study Area and Crop Management

The study was carried out in the Colombian Andean region, on the high plateau in the department of Cundinamarca, in the municipality of Subachoque (4.888668° N, 74.18668° W, ~2609 m above sea level) (Figure 1, left). The EC footprint is also indicated, enclosed in the red curve in Figure 1, right.

This region has soils that are apt for agriculture, and there is water availability for crop irrigation, supporting the predominant economic activity, agriculture, given its proximity to the capital of the country, Bogota D.C., which is home to 7,834,167 habitants (2021), with a continuous high demand for food, such as vegetables. The soil is part of the Toscana soil unit [35], represented by the taxonomic soil name Typic Hapludands. The landscape is a fluvial–lacustrine plain, with terraces at lower elevations and flat slopes (0–1%). The soils are derived from volcanic ash on coarse fluvial materials. From the climatic point of view, its average multiannual precipitation ranges between 750 and 1000 mm (with a bimodal distribution), average annual air temperature ranges between 12 and 14 °C, multiannual relative humidity is between 80 and 85%, and sunlight ranges from 1400 to 1500 h [36]. The carrot production system was evaluated in a commercial plot of 3.11 ha, under fixed sprinkler irrigation, using three seed varieties, as follows: 50% Olimpo de Saenz Fetty, 25% Cardiff, and 25% Miraflores, with a sowing density of 480 g of sexual seed per hectare at 8 cm (about 3.15 in) distance between plants and 50 cm (about 1.64 ft) between furrows. The planting date was 10 August 2021, and the harvest date was from 25 December 2021 to 22 January 2022, a period in which there was an ENSO (El Niño–Southern Oscillation) climatic phenomenon in its cold phase, known as La Niña, where historically there is an increase in precipitation, a factor that largely regulates the production cycle of vegetables in the region. The crop was grown following the common cultural management that includes foliar and soil fertilization, weed control, insecticide, and fungicide treatment.

2.2. Microclimate and EC System Measurements

The EC Flux tower system was deployed at a farmer’s farm, on a plot that has been using conventional tillage practices. The EC technique continuously recorded net carbon ecosystem exchange (NEE) and microclimatic variables. The equipment was installed on 16 August 2021 (6 days post-planting, DPP), and the measurements were taken until the end of the harvest on 20 January 2022 (157 DPP). The EC tower consists of an IRGASON with a gas analyzer of an open path (EC 150, Campbell Scientific, Inc., Logan, UT, USA) and a 3D sonic anemometer (CSAT3A, Campbell Scientific, Inc., Logan, UT, USA). The IGARSON and 3D anemometer are operated using a separate electronic module (EC100, Campbell Scientific, Inc., Logan, UT, USA). The EC tower used here and the configuration of each sensor is described in [35,36].

NEE Partitioning

Starting from the premise that NEE = GPP + Rd [37] and according to [38,39], the non-linear Mitscherlich light response function represents the net ecosystem exchange (NEE) as a function of the photosynthetic photon flux density (PPFD). This function was used to partition the diurnal NEE (solar global radiation > 1 W m⁻²) into ecosystem respiration (R_eco) and GPP [38,39], as follows:

$$NEE=-\left(A_{gmax}+R_d\right)\times \left(1-{\exp }^{\left(\frac{-\Phi *I_{inc}}{A_{gmax}+R_d}\right)}\right)+R_d$$

(1)

where $R_d$ is the respiration term [µ_mol (CO₂) m⁻² s⁻¹], $A_{gmax}$ is the CO₂ uptake at light saturation [µ_mol (CO₂) m⁻² s⁻¹], $\Phi$ is the quantum yield efficiency [µ_mol (CO₂) µ_mol (photon) ⁻¹] corresponding to the initial slope of the photosynthetic light response curve, and $I_{inc}$ is the incident PPFD [µ_mol (photon) m⁻² s⁻¹]. For each day, the set of parameters were calculated using non-linear regression and a subset of NEE and PPFD data, within a moving window of 14 days, centered on each day. For each daily diurnal half-hour, GPP was estimated by subtracting R_d from the non-linear Mitscherlich light response function, as follows:

$$\mathrm{GPP}=-(A_{gmax}+Rd)\times (1-{\exp }^{\left(\frac{-\Phi *I_{inc}}{A_{gmax}+R_d}\right)})$$

(2)

Here, R_eco was calculated by subtracting the GPP estimated from the measured NEE. Quality control, gap-filling, data postprocessing, energy balance closure, uncertainty analysis, and statistical analysis methods were applied, following the procedures described in [35]. The daily GPP reported here is measured in g C m⁻² day⁻¹.

2.3. Multispectral Imagery

Multispectral imagery of the farm was acquired with a Sequoia Parrot multispectral camera mounted on the gimbal of an aerial platform DJI Phantom 3 Advanced quadcopter drone, flying at a height of 30 m above the field. The camera acquires four spectral bands in the green (530–570 nm), red (640–680 nm), red-edge (730–740 nm), and NIR (770–810 nm) bands, plus an RGB image to help the registration process. The images have a resolution of 1.2 megapixels and a radiometric resolution of 16 bits. This camera has a light sensor that records solar irradiance data for each spectral band and which was used for sunshine calibration. A campaign of image acquisition was conducted on the field on 11 dates in 2021, namely 7 September 2021, 23 September 2021, 30 September 2021, 8 October 2021, 15 October 2021, 22 October 2021, 1 November 2021, 10 November 2021, 15 November 2021, 23 December 2021, and 27 December 2021, starting with the emergence of the crop until its senescence. Due to atmospheric conditions some flights had to be canceled, resulting in a non-uniform distribution of the flight dates. Besides the images taken with the UAV over the farm, an extra image of a Spectralon reference panel placed on the ground was acquired with the multispectral camera to calibrate radiance images to reflectance images. No atmospheric correction is necessary, since UAV data is far less influenced by atmospheric effects. Before reflectance calibration, a radiometric and geometric correction of the four spectral bands was performed in Python using the specialized code for the Sequoia camera found on GitHub [40]. After radiometric and geometric corrections, the MicaSense software also allowed us to perform sunshine calibration by dividing the radiance of the soil by the radiance of the sun. The conversion to reflectance is performed using the following equation:

$$\rho \left(x,y,\lambda \right)=\frac{I\left(x,y,\lambda \right){\rho }_S\left(\lambda \right)}{Is\left(\lambda \right)}$$

(3)

where $\rho \left(x,y,\lambda \right)$ is the reflectance image at pixel coordinates $x,y$ and waveband $\lambda$, $I\left(x,y,\lambda \right)$ is the raw intensity image at pixel coordinates $x,y,$ and waveband $\lambda$, ${\rho }_S\left(\lambda \right)$ is the known reflectance of the Spectralon panel at $\lambda$ wavelength (0.99 at visible and NIR ranges), and $Is\left(\lambda \right)$ is the mean intensity of the Spectralon panel at waveband $\lambda$ [41].

Once the images of the drone are calibrated, an orthophoto mosaic is constructed for each band using six ground control points in AgiSoft software. The four orthophoto mosaics corresponding to each spectral band are reprojected to have the same size, with a pixel size of 3.11 × 3.11 cm (about 1.22 in). Using the shapefile delimiting the footprint of the EC, a clip of each orthophoto mosaic is made, containing the region of the crop covered by the EC footprint, which provides a first mask to delimit the region to compute the VIs. A second mask is made to differentiate plants from the soil, so that the VIs are computed on vegetation pixels only. We use the NDVI to differentiate plants from the soil by setting a threshold of NDVI > 0.3 to detect plants, based on [42] and prior experience. The intersection of the two masks is used to compute the VIs, within the EC footprint.

2.4. Spectral Indexes

A total of 24 VIs were studied to determine their capacity to predict the daily GPP estimated the EC method. We implemented all of the indexes studied in [43] using the four bands available, namely NDVI, SAVI also referenced in [1], the modified soil adjusted vegetation indexes (MSAVI, MSAVI2), transformed normalized difference vegetation index (TSAVI), MTCI, also refenced in [39], the infrared percentage vegetation index (IPVI), transformed normalized difference vegetation index (TNDVI), green normalized difference vegetation index (GNDVI), global environmental monitoring index (GEMI), pigment specific simple ratio (PSSR), ratio vegetation index (RVI), MCARI, also referenced in [20], the difference vegetation index (DVI), and WDVI, as referenced in [1]. We could not implement EVI, because it requires a blue band. However, we implemented the enhanced vegetation index 2 (EVI2) [44] that closely follows EVI using only the red and NIR bands. In addition, we also implemented CIgreen [17,18], CIred-edge, and red-edge NDVI (RENDVI) indicated in [20]. We also included the normalized difference red edge index (NDRE) found in [45] to add more red-edge-based indexes, WDRVI [17], and WDRVIgreen [18]. Finally, we included two novel VIs that have been found to correlate with GPP, namely the green–red vegetation index (GRVI) [29] and near-infrared reflectance of vegetation (NIRv) [19]. We did not find good results for TSAVI using the parameters indicated in [43], which coincides with other papers using the TSAVI; instead, we used the predetermined values indicated in [46] that use a stronger X parameter to suppress noise, which produced better results.

After studying the graphs between the daily GPP and Vis, we found that linear regression was a good model to represent the relationship between the VIs and the GPP on several indexes (see Figure 2). In most VIs, the relationship GPP–VIs does not follow any discernible relationship or correlation, and we employ linear regression to be consistent with all the other indexes (Figure 2). The root mean square error (RMSE) and the R² score were computed to estimate the performance of the linear regression on each index, as follows:

$$RMSE=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^N{\left(x_i-{\widehat{x}}_i\right)}^2}{N}}$$

(4)

$$R^2=1-\frac{{{\displaystyle \sum }}_{i=1}^N{\left(x_i-{\widehat{x}}_i\right)}^2}{{{\displaystyle \sum }}_{i=1}^N{\left(x_i-\overline{x}\right)}^2}$$

(5)

where, N is the number of points, $x_i$ the observed value, and ${\widehat{x}}_i$ the value predicted by the VI.

The implementation of all the methods reported here was carried out in Python 3.8. Linear regression was performed using the ordinary least square (OLS) model in the statsmodel module. The RMSE and R² statistics were computed using the mean_squared_error and r2_score functions from the sklearn.metrics module. We perform a linear regression of the GPP vs. VIs using all the data available to screen out the best indexes using the RMSE and R² statistics, since we are not interested in the model parameters of all models, only in the models that perform better. To predict the performance of the best models on unseen data, we use leave-one-out cross-validation (LOOC), since the dataset is small (11 dates of image acquisition), and it has been shown that LOOC is an almost unbiased estimator of the true model generalization performance [47]. We select the first data sample as the test sample consisting of a VI value on the first date and the corresponding GPP value from the EC. The remaining 10 samples are used to train the first linear model, which is used to compute the error on the test sample. Then, the second data sample is selected as the test sample, and the remaining 10 samples are used to train a second linear model, which is used to compute the error on the second test sample. The process is repeated until all samples are used as test samples. In the end, we end up with 11 linear models and the estimated error on each of the 11 samples. This error is an estimate of the generalization error of the linear model. The error on each test sample corresponds to $x_i-{\widehat{x}}_i$ on Equations (4) and (5). Hence, we can compute the cross-validation RMSE and R² statistic. In addition, we also have 11 linear models with parameters varying from each model to another, since they are trained on a slightly different data set. Hence, we can choose the average of the two model parameters, namely intercept and slope, as a good estimate of the linear model parameters. An estimate of their standard error (BSE) can be computed from the dispersion provided by the cross-validation. In addition, we can compute the p-value of the mean coefficients using a t-test, since the dataset is small. The t-values of this test consist of dividing the mean coefficient by its BSE. Then, the p-value can be obtained from the t-student distribution with 9 degrees of freedom (11 points minus 2 model parameters). This p-value can be obtained in Python from the scipy.stats module.

3. Results

Figure 2 shows the GPP–VI samples, the regression line, and the RMSE and R² statistics. As can be seen from this figure, some VIs follow a linear trend more closely (SAVI, MSAVI, MSAVI2, TSAVI, GEMI, DVI, WDVI, and EVI2), others show moderate correlation (NDVI, IPVI, TNDVI, RVI, and WDRVI), while the remaining indexes show no correlation with the daily GPP.

Figure 3 shows more compactly the RMSE and R² statistic of all the VIs. This figure allows us to extract the top seven VIs that correlate well with the GPP; they are SAVI, MSAVI, MSAVI2, DVI, WDVI, EVI2, and NIRv. It is a surprise that the simple DVI provides good regression results. Furthermore, it is noteworthy that all VIs based on the red-edge (NDRE, PSSR, RENDVI, CIred-edge) and green (GNDVI, Cigreen, GRVI) do not correlate with the GPP. This could mean that the red and NIR bands suffice to obtain a good estimation of the GPP in carrot crops.

We selected the best models found to predict the daily GPP from VIs, which can be used in further studies of the GPP in carrot crops, under different environmental conditions and locations. The best VI was NIRv, with the lowest RMSE and highest R² score (0.8465), followed by MSAVI (0.8323), MSAVI2 (0.8317), SAVI (0.8257), EVI2 (0.8256), WDVI (0.8207), and DVI (0.8154).

The results of the cross-validation are indicated in

Table 1. Models of the selected VIs to predict the daily GPP with their cross-validation statistics.

Vegetation Index	Model	RMSE	R2	p-Value
NIRv	GPP = −0.28(0.16) + 24.12(0.99) × NIRv	1.6589	0.7823	0.06	7.8 × 10−10
MSAVI	GPP = −1.37(0.21) + 16.96(0.68) × MSAVI	1.7780	0.7667	4.7 × 10−5	6.5 × 10−10
MSAVI2	GPP = −1.28(0.20) + 15.71(0.64) × MSAVI2	1.7898	0.7652	7.2 × 10−5	7.6 × 10−10
SAVI	GPP = −2.09(0.24) + 17.55(0.72) × SAVI	1.8446	0.7580	6.3 × 10−6	7.7 × 10−10
EVI2	GPP = −1.45(0.22) + 15.45(0.67) × EVI2	1.8730	0.7543	4.6 × 10−5	1.3 × 10−9
WDVI	GPP = −0.56(0.16) + 20.38(0.67) × WDVI	1.8091	0.7626	3 × 10−3	1.1 × 10−10
DVI	GPP = −0.58(0.15) + 20.06(0.64) × DVI	1.8468	0.7577	2.4 × 10−3	7.9 × 10−11

. The model indicates the mean intercept and the mean slope, and the BSE of these parameters is reported in parenthesis. The cross-validated RMSE and R² of the model are also indicated in the table. The p-value column contains the p-value for the intersect on the left and the p-value for the slope on the right column.

With these seven indexes, we want to estimate the daily GPP and compare it with the EC-estimated daily GPP. Since we only have 11 dates and the image acquisition was from September 7 to December 27, i.e., 112 days, we need to interpolate either the VIs and use the models to estimate the daily GPP or to interpolate the GPP estimated on the 11 dates to the remaining dates. Here, we choose to interpolate the GPP directly, since simple linear interpolation of the Vis’ time series leads to a linear time series of the GPP, which is incorrect. Figure 4 shows the estimated GPP as a function of the day of the year, using the NIRv index, which is the best (Table 1). This figure suggests that we can linearly interpolate the missing dates. Notice that this is not a model with parameters fitting the 11 points; rather, it is a simple interpolation between each pair of consecutive points. No cross-validation can be performed, since a LOOC will try to fit one single point from the remaining points and we do not want to estimate the points that we already have with a cross-validated model; instead, we want to estimate the daily GPP from the points we have. We could build a polynomial model to fit the points, but since there is a strong variability (Figure 4) the polynomial will be of a high order (high complexity), and we cannot report a polynomial model with its parameters cross-validated using just 11 points. It would be wrong to report a model that predicts GPP on any day of the crop growth cycle, since it would mean we no longer need VIs, just the day of the year. What we are proposing here is not a model to estimate the daily GPP, but rather an easy (model-free) way to interpolate the missing dates, with no training required. This simple strategy can be used on other datasets, which will most likely have different dates and probably a different behavior than the one observed here and, hence, we cannot evaluate this strategy as a model with parameters.

Figure 5 shows the temporal trajectory of the estimated GPP using simple interpolation for all the indexes selected, with their respective RMSE and R² score (No BSE can be computed).

The accumulated GPP is an important measure that is used to predict yield and efficiency in the use of resources (water, carbon, and nutrients). The accumulated GPP-EC is 580.4 g C m⁻², and the estimated accumulated GPP using NIRv is 520.4 g C m⁻² (10.3% error); using MSAVI it is 516.7 g C m⁻² (11.0% error); using MSAVI2 it is 515.3 g C m⁻² (11.2% error); using SAVI it is 514.9 g C m⁻² (11.3% error); using EVI2 it is 513.6 g C m⁻² (11.5% error); using WDVI it is 521.5 g C m⁻² (10.1% error); and using DVI it is 522.1 g C m⁻² (10.0% error). Surprisingly, the estimated accumulated GPP is a bit better using WDVI and DVI even though the other indexes outperform them. As indicated in [10], using a more complex LUE model with an R² score of 0.96, the error in the accumulated GPP is around 6%, so the 10% error found here can be expected from a much simpler model.

4. Discussion

The moderate performance of the NDVI index to directly predict GPP is a surprise, since it has been shown to perform well using satellite data [10,13,16]. However, in [1,18] the direct NDVI-GPP correlation is moderate. A Landsat-7 study [17] shows good performance of NDVI, but also indicates that the NDVI exhibits limitations in GPP estimation at moderate to high vegetation densities. In [1], a study in maize using MODIS imagery shows that VIs (NDVI, SAVI, EVI, WDVI) are moderately linearly related to GPP, with NDVI being the worst. This study also shows that a good correlation between the VIs and the GPP can be obtained using a non-linear relationship including the PAR, as follows: VI × VI × PAR. Related research [20] uses VI × PAR with moderated correlations for remote sensing images. In [10], another study in maize using Landsat-8 shows better results using the VPM. The EVI has been amply used to estimate GPP with good correlation [10,14,17,18,21] as well as the EVI2 [28], coinciding with our good results here with EVI2. A moderate correlation of EVI with GPP has been shown in [1], and good-to-moderate performance was shown in [18]. We obtain a good correlation between SAVI and GPP, and similar to the results indicated in [1]. We also obtain a good correlation between MSAVI and MSAVI2 to GPP, with MSAVI being the best. The MSAVI indexes have not been reported in the literature as showing a good correlation with GPP, so they should be included in further studies. The simple DVI and WDVI showed good correlation with GPP on the 11 spectral samples and, surprisingly, were a bit better than the best index NIRv for estimating the accumulated GPP. The novel NIRv has been shown to have a very good correlation with GPP in [19,28], and was the best VI at the 11 temporal samples taken. It is interesting to note that the NIRv is the product of the scene NIR and the NDVI, and that this simple correction of the NDVI overcomes the limitations of the NDVI, agreeing with our results. Contrary to the results indicated in [17,18,20], the CIgreen, the CIred-edge, the RENDVI, the MCARI, the MTCI, and the WDRVIgreen show no correlation at all with GPP. As indicated in [1,2,10,12,16,18,20,21,22], empirically-based VI models using additional information to the VIs provide, in general, better correlations between VIs and GPP than VIs alone, and LUE-based models perform even better than empirically-based VIs. Our results are good compared with the general results using satellite imagery and simple vegetation indexes. However, they can be improved using more complex models that can obtain R² scores well above 0.9 to better predict the daily GPP during a crop’s lifecycle.

These results are an important first step to being able to routinely monitor carbon uptake on carrot crops on regions that are highly covered by clouds. Other crops which are highly cultivated in the area, such as potatoes, could also be studied to develop GPP–VIs linear models, as we accomplished here, to monitor carbon uptake using high-resolution, cloud-free multispectral imagery acquired with UAVs.

5. Conclusions

This work shows that it is possible to estimate the GPP of a carrot crop using high-spatial resolution UAV-based multispectral imagery and simple linear regression models, which can be used to monitor GPP on small or medium size farms that cannot be monitored using satellite imagery due to coarse resolution and high cloud cover present in humid areas, as is the case in Colombia. The integration of cloud-free high-resolution UAV imagery and the EC method can be used to obtain accurate models to estimate carbon uptake from VIs in several agroecosystems. In the context of precision agriculture, estimating GPP from VIs computed at high spatial resolution while avoiding mixing soil and vegetation is an important step to estimate crop yield, and to ensure the efficient use of resources (water, carbon, and nutrients).