A NOVEL METHOD FOR ESTIMATION OF STRUCTURAL CHANGES IN POTATO CROPS FROM UAV-BASED DIGITAL SURFACE MODELS

Climate change is arising challenges in potato crop production. More robust clones and cultivars adapted to climate change and resistant to pests and diseases, need to be developed. Several phenotyping approaches are being explored for a quick and affordable characterization of structural changes of crop varieties in real-world conditions. However, conventional methods present several problems. Field work in remote areas is time-consuming or difficult to carry out because of physical barriers. Satellite images lack the temporal resolution needed to monitor small changes in plant structure. Unmanned aerial vehicles (UAVs) provide highresolution images which can give a better understanding of vegetation changes. Recent UAV-based studies have demonstrated the potential of the Normalized Difference Elevation Index (NDEI) to map centimeter-level changes in topography in different environments. This paper explores the potential of NDEI for identifying structural changes in experimental potato crops. In particular, it proposes an automatic method for estimation of volume changes in potato crop canopy. NDEI values were obtained by multitemporal analysis of digital elevation models generated from UAV images using structure-from-motion (SfM) and multivision-stereo (MVS) techniques. The main components of the proposed method are: i) Image data capture from UAVs, ii) Data pre-processing, iii) NDEI calculation, and iv) Data processing and analysis. The results suggest that changes in surface elevation values, normalized using the NDEI index, are useful for rapid assessment of volume change in potato crop canopy. The NDEI reveals as an affordable technique to field measurement of phenotypic traits in potato breeding programs.


INTRODUCTION
Potato (Solanum tuberosum) is the fourth most important crop in the world after maize, rice and wheat (Namugga et al., 2018). In Colombia, the largest commercial potato industry in the Andean region, it is estimated that around 400,000 families are involved in this industry, from harvesting to commercialization (Dane, 2017). Colombian potato crops are grown between 1,800 and 2,500 meters above sea level depending on varieties and species. Colombian potato production is concentrated in three departments: Cundinamarca (45.54%), Boyaca (25.08%), and Nariño (19.28%) (Dane, 2017).
In Colombia, as well as in other countries, agricultural production is facing pressure to meet the societal challenge of jointly ensuring food security, adaptation to climate change impacts and mitigation of greenhouse gas emissions (Aalborg-University-Copenhagen, 2017). This challenge calls for prompt action from producers and government agencies to make improvements to current agricultural practices. Moreover, there is a need to breed robust plants, adapted to climate change and resistant to pests and diseases. As plants are intrinsically related to the surrounding environment, observed phenotypes are a direct product of this interaction. Therefore, the ability to study and quantify phenotypes in the field, under real-world conditions, is essential to the basic understanding and improvement of structural traits (Pauli et al., 2016).
In order to develop successful potato breeding programs, it is critical to identify low-cost technologies that help to accurately measure phenotypic traits of a variety of clones and cultivars. * vdangulom@correo.udistrital.edu.co Unmanned aerial vehicles (UAV) is a low-cost technology which may carry remote sensors to support field phenotyping studies. One of the main advantages of using UAV-based data in agricultural research is the ability to perform high-resolution temporal sampling of crops and to generate high spatial resolution digital surface models (Adão et al., 2017, Deng et al., 2018, Koh et al., 2019 The analysis of structural vegetation changes with remote sensing includes different techniques, the main methods being: (i) mapping of changes on canopy reflectance; and (ii) mapping of changes on 2D morphological traits (Pauli et al., 2016). Canopy reflectance and 2D morphological traits can provide insight into overall plant health as well as specific physiological processes. To date, most field-based phenotyping and remote sensing applications have focused on using multispectral vegetation indices to infer overall plant status (Pauli et al., 2016), with normalized difference vegetation index being the most well known index.
Although 2D indices and 2D morphological traits generally can be informative, they lack the ability to provide vertical information needed to understand tridimensional processes. In contrast, approaches such as structure from motion (SfM) and multistereo vision (MSV) provide point cloud data which may be used to obtain three-dimensional information of crop morphological traits, both for crop canopy at plot level and for individual plants' height, shape and volume (Paulus et al., 2014) Recent studies have demonstrated the potential of the Normalized Difference Elevation Index (NDEI), based on UAV images, to identify changes in land surface elevation in different environments (Morales et al., 2019, Lizarazo et al., 2017, Zhang et al., 2019. However, the potential of NDEI for deriving 3D crop structural information in agricultural studies has not been yet evaluated. The objective of this paper is to develop a novel technique that allows estimating structural changes in potato crops from digital surface models. The representation and analysis of changes in potato crop canopy is done using the NDEI index. The method proposed is tested for detection of volume changes in potato crop canopy in an experimental field in Cundinamarca (Colombia).

DATA AND METHODS
The method proposed in this paper uses the NDEI index obtained from UAV-based multispectral images to map structural changes in a potato crop field phenotyping experiment. The NDEI was applied for identification of individual plants which were intentionally removed from the crop to calculate biomass as well as to estimate changes in crop canopy volume due to such removal. The procedure for analyzing structural changes is: (i) data capture from UAV; (ii) data pre-processing using SfM-MVS techniques; (iii) NDEI calculation; and (iv) data processing and analysis.

Study area
The study area is an experimental potato crop of approximately 2 ha, located in an agricultural zone in Mosquera, Cundinamarca (Colombia). The study area is 14 km west to the city of Bogotá with average altitude of 2560 meters above sea level. The study area is part of the Universidad Nacional's Marengo agricultural research centre which covers approximately 97.4 ha (See Figure 1).
The UAV flight mission was designed to handling a 80 percent forward overlap and a 60 percent lateral overlap, and a cross fligh at 40 meters altitude above ground level. The UAV mission comprised 14 flight lines to acquire 500 images and had a flight time of 8 minutes at a speed of 5 meters per second. Multispectral images were acquired from 11:00 am to 1:00 pm, Colombia local time . This ensured the necessary lighting conditions to minimize the effects of shadows on the crop. This schedule was chosen to avoid the strong winds characteristic of the area during the afternoon hours. (See Figure 2 Four flights were made, two per day, on 28th December 2019 and on 4th January 2020. The photogrammetric flight routes were programmed using the ground control station software. Captured images have a spatial resolution of 2.6 cm/pixel. In order to obtain high positional accuracy digital elevation models, a total of 7 markers were measured using high-precision surveying instruments to be used as Ground Control Points (GCP) and Control Points (CP). Figure 3 shows the design of the GCP marker, a fiberglass device built to cover a 45cm side square. The GCP design was made in the form of black and orange squares. The size of the GCP was selected to make sure it covers at least 30 pixels in a UAV image acquired with a Micasense RedEdge sensor at a height of 100 meters above ground.

Data Pre-processing
The UAV-based captured images were radiometrically corrected with a Spectralon calibration panel to obtain pixel values between 0 and 1 representing ground suface reflectance and geometrically corrected through a photogrammetric adjustment process which involved the integration of GCPs and CPs. The photogrammetric adjustment of the UAV images was carried out using the Agisoft Photoscan software. The adjustment finds homologous points by analyzing all the images captured with the UAV. The algorithm uses the Scale Invariant Feature Transform (SIFT) for extracting distinctive invariant features from images to perform reliable matching between different views of the scene (Lowe, 2004) (See Figure 5).
Georreferencing was achieved by using GCPs during the photogrammetric processing, CP were used to assess model accuracy for flights 1 and 2 (See Figure 4). Several issues during data capture affected the usability of two markers on 4th January 2020, leading to avoid CP in flights 3 and 4. Table 1 shows number of GCPs and CPs used for each flight.  Image pre-processing generates two outputs: (i) an orthorectified mosaic, also known as orthophoto; and (ii) a digital elevation model (DEM, See Figure 6). These two outputs were used in the data processing stage to obtain structural changes in potato crops from the normalized elevation difference index (NDEI).

NDEI: Normalized Difference Elevation Index
The Normalized Difference Elevation Index (NDEI) is a metric that can be used to identify and quantify changes in digital elevation models from UAV images. Two digital elevation models obtained in different time periods, are needed to calculate the NDEI index and thus to identify elevation changes occurred : where DEM = Digital elevation model's value DEMt1 = Value of the early DEM DEMt2 = Value of the later DEM The NDEI index between two time periods represent a decrease in surface elevation when values are lower than 0, and an increase in surface elevation when values are greater than 0. In crop fields, surface elevation changes are usually caused by canopy growth; however, they may also correspond to vegetation removal. The latter case is exploited in the method proposed here.

Data processing and analysis
In this stage, the analysis of NDEI values was conducted. It starts with the digital elevation models (DEMs) calculated on the dates of UAV data collection. Each consecutive pair of DEMs were used to calculate NDEI for the evaluation of structural changes in the potato crop during the period.
To identify zones of structural change from the NDEI indices, the Python programming language and the skimage library were used ( Van der Walt et al., 2014). Several algorithms provided by this library were used to process the NDEI indices, as described in the next paragraph (Lizarazo et al., 2017). • Threshold method: It computes the optimal threshold in the histogram of a raster object, maximizing the variance between two kinds of pixels, which are separated by the threshold. This method is applied to the raster object in order to have a value that allows separation of the objects from the background obtaining two values (binary object). • Closing method: This morphological method of an image is defined as a dilation followed by an erosion. It allows to close objects that have holes inside itself. • ClearBorder method: This method allows to eliminate objects that are incomplete in the edges of the image, these elements must be eliminated since their contents does not meet the requirements and characteristics to be identified. • Label method: This method allows to join neighboring pixels as a single object in an array. In this function we create an array of regions that are classified as plants. • The volume calculation is done by averaging the elevation values of each pixel, in meters, that belongs to the object detected in the digital elevation model; then, the resulting value is multiplied by the size X and Y of pixel; the output of the algorithm is a value in m 3 . Figure 7 shows the result of the algorithm implemented to identify 3D changes, i.e. to detect crop zones where abrupt changes in height and surface happened in a short period of time due to removal of individual potato plants. Each detected change zone is labeled with a number and a red box, this allows easy identification of the change. Figure 8 shows a zoom on one of the change zones; the pixel values correspond to the NDEI values in the change zones.
On December 28th, field work was carried out to remove four potato plants and perform a destructive biomass sampling of tubers and leaves. Figure 10, shows the orthophoto corresponding to the flight carried out on December 28 at 1:25PM immediately after removing the plants. The plant removal sites, that is, the areas of surface elevation change in crop canopy, correspond to the red colour boxes, are identified with IDs 1, 2, 3, and 4. Figure 11 shows the DEMs obtained for the 4 flights: a) December 28, 2019 at 12:10PM; b) December 28, 2019 at 1:25PM; c) January 4, 2020 at 12:01PM; d) January 4, 2020 at 12:40PM. It can be seen several differences in crop canopy surface elevation between the time periods. However, an accurate visual identification of zones of elevation change from the DEMs is a complicated task.
For each flight, the root mean square error (RMSE) of the GCPs was calculated and used to evaluate the geometric quality of the model. Table 2, Table 4, Table 6, and Table 7 Table 3 summarizes accuracy of DEM for the flight 1 on December 28th, 2019 at 12:10PM (GMT-5), markers 5 and 6 were used as CPs. It can be seen that the highest RMSE value is 3.52 cm that corresponds to the Z axis, the lowest RMSE value is 0.35 cm and corresponds to the X axis. Total accuracy was 3.71 cm. RMSE of the DEM for the flight 2 on December 28th, 2019 at 1:25PM (GMT-5) can be seen in the Table 5, marker 6 was used as CP, the highest RMSE value is -6.88 cm and corresponds to the Z axis as well as flight 1. The lowest RMSE value is -1.28 cm and corresponds to the X axis. Total accuracy for this flight was 7.30 cm.
Problems with two capture devices limited the number of markers available for the photogrammetric processing on the January 4th flights, preventing from using CPs on these flights as the minimum number of markers to achieve a good quality model were 5. It can be seen in Table 6 and Table 7 that the RMSE values are unusually high (RMSE values are shown in mm). The photogrammetric processing of these flights was carried out    Figure 12 shows the NDEI index calculated: (a) NDEI index calculated for the period: January-December or flight 4 -1. Red colour indicates areas where elevation decrease is present (i.e. there is a loss in crop canopy volume) ; green colour represent areas where elevation increase is minimal or zero. (b) Flights performed on December 28, 2019 (Flight 1 and 2); (c) Flights performed on January 4, 2020 (Flight 3 and 4). Red colour indicates areas where either a negative difference between the two periods exists, or a loss of volume happens. Green colour indicates areas where there is no surface elevation change or minimal change.
The NDEI index of December 28 shows four zones of volume loss on the potato crop canopy. These changes correspond to the destructive removal of potato plants to estimate individual biomass. This change in crop canopy can be identified in the NDEI index map because the removal task was done in the middle of the two flights, that is, at some time between 12:10PM and 1:25PM. The NDEI index map also shows small changes in the crop canopy surface elevation due to the pedestrian traffic of two researchers who removed the potato plants. None of the above changes can be easily identified neither in the orthophotomosaics nor the digital surface models. The NDEI index for January 4 shows no significant changes, due that no field work was conducted on this date that would represent crop surface elevation canopy changes between the two flights performed on this day.
Algorithm 1 is applied to automatically identify the changes from the NDEI index. The implemented algorithm takes out a matrix of pixels and outputs tagged change zones. Figure 13 shows the result of the automatic identification of the change zones from the NDEI index: (a) It shows the zones of change calculated from the NDEI of Flights 4 and 1 which is an index calculated between two dates January 4 and December 28; (b) It shows the zones of change calculated from the NDEI of Flights 2-1 which is an index calculated between two flights of the same date December 28; and (c) It shows the zones of change calculated from the NDEI of Flights 4 and 3 which is an index calculated between two flights of the same date January 4.
The implemented algorithm is capable of estimating the volume of the identified change zones. Table 8 summarizes data for the objects automatically detected from each NDEI index evaluated,. It also shows total volume and average volume (m 3 that changed between the two time periods under study. In the potato crop plot, four plants were removed on December 28, 2019. Table 9 shows the volume of change in each NDEI index in the four zones corresponding to the plants. It is important to note that: (i) in the NDEI index of January 4 between the two flights made on this day, the estimated change is close to zero; and (ii) the detection algorithm does not find a significant change on this day; consequently, the estimated change values are also 0 for the four zones.   Figure 14, shows the automatically identified change zones where individual plants were removed during the field work. It can be seen that between the two dates there is a decrease in the NDEI values in the four zones where the plants were removed for field work. This decrease in the NDEI value is associated with the volume values calculated by the algorithm. Corresponding to a change that occurred due to movement or removal of elements in the study zone over a period of eight days.
These results demonstrate that the proposed method can be used to estimate structural changes in crop canopy and volume. However, it should be noted that the potato crop yield potential is dependent on fast and homogeneous canopy growth, early attainment of maximum leaf area for maximum light intercep-   Table 8. Summary of volume changes obtained by applying the implemented algortithm to the NDEI indexes tion, and completing growth within the available growing season. Thus, the proposed method may be further used in future field-based phenotyping studies to estimate both 2D and 3D morphometric variables.

CONCLUSIONS
The NDEI index was applied in a potato crop experiment to detect 3D structural canopy changes in different time periods. We analyzed two time periods with a difference of 8 days, in order to extract relevant information related to structural changes in potato crop canopy surface elevation and volume. Qualitative evaluation of the proposed methods, in terms of identification of surface elevation change zones, suggest a good capacity for  It should be noted that, for successful application of the proposed method, the accuracy of the digital elevation models should be centimeter-level because even small changes in elevation values may represent significant changes in the output NDEI indices. This means that, in order to isolate actual changes in elevation from bogus changes due to positional error, it is very convenient to use high precision surveying instruments to de-termine GCPs position.
The proposed method, implemented in Python, aimed to determine surface elevation changes without manual intervention. It applies computer vision techniques to identify zones of change.
The accuracy of the proposed method seems to be very competitive compared to traditional manual field techniques which are time-consuming and expensive.
Although the proposed approach is still experimental, the results of this study suggest that it is feasible to estimate the area and volume of crop canopy changes from images captured by UAV cameras. The proposed approach may help farmers and breeders to have low-cost indicators of crop canopy change who may be useful for agricultural decision making.
UAV-based photogrammetric remote sensing is an option for estimating structural and volume changes in potato crops and, thus, understanding canopy crop development and growth. This paper's results confirm previous studies that have shown that multi-temporal remote sensing allows affordable crop monitoring for smallholder farmers. Systematic use of methods like the one presented here allows producing high spatial digital surface models to derive phenotypic traits relevant for crop adaptation to changing climatic conditions.