The Outlining of Agricultural Plots Based on Spatiotemporal Consensus Segmentation

The outlining of agricultural land is an important task for obtaining primary information used to create agricultural policies, estimate subsidies and agricultural insurance, and update agricultural geographical databases, among others. Most of the automatic and semi-automatic methods used for outlining agricultural plots using remotely sensed imagery are based on image segmentation. However, these approaches are usually sensitive to intra-plot variability and depend on the selection of the correct parameters, resulting in a poor performance due to the variability in the shape, size, and texture of the agricultural landscapes. In this work, a new methodology based on consensus image segmentation for outlining agricultural plots is presented. The proposed methodology combines segmentation at different scales—carried out using a superpixel (SP) method—and different dates from the same growing season to obtain a single segmentation of the agricultural plots. A visual and numerical comparison of the results provided by the proposed methodology with field-based data (ground truth) shows that the use of segmentation consensus is promising for outlining agricultural plots in a semi-supervised manner.


Introduction
The world's population has tripled over the last 100 years and is still growing dramatically while resources have remained the same, causing changes in the outlook of food supply.According to the Food and Agriculture Organization (FAO), global food production will need to grow by 70% in order to satisfy the food and feed demand of a population of 9 billion people by 2050.The agriculture sector faces a critical overall challenge: to ensure access to safe, healthy, and nutritious food while using natural resources more sustainably and making an effective contribution to climate change adaptation and mitigation [1].
This challenge implies a greater pressure than ever before on productive land.Accurate and up-to-date information about agricultural land, such as its status, acreage, ownership, and the type of crops, allows stakeholders to establish effective agricultural policies (e.g., for reducing greenhouse gas emissions, regulating water rights, and estimating subsidies and agriculture insurances) [2,3] and update agricultural geographical databases [4], among other important tasks.In order to have up-to-date information on agricultural land, it is essential that the outlines of the plots are correct and can be quickly updated.
For the last 40 years, the outlining of plots in agricultural lands has been addressed through different initiatives around the world.In the United States, the National Research Council has published two reports that examine the situation of the land parcel data in the U.S. and provide a series of recommendations that would foster a national data system for storing plot information [5,6].On the other hand, the European Union has promoted the use of a common Land Parcel Identification System (LPIS, https://ec.europa.eu/jrc/en/scientific-tool/lpis-quality-assessment)among its members in order to maintain a record of the activities of farmers on their lands [7].The success of all these initiatives depends enormously on a precise outlining of the parcels.
Traditionally, the outlining of the parcels has been carried out manually using photointerpretation or field campaigns and documented by surveying sketches and textual documentation, all of which is very expensive both in time and financial resources.
The application of Information Technology tools, such as remote sensing and geographical information systems, improves efficiency in the agricultural sector, enabling planning and decision making based on the spatially and temporally distributed data provided by these tools [8,9].The new generation of optical remote sensors placed on aircraft, satellite platforms, and drones offers accessible and useful data of very high resolution for monitoring agricultural fields at a plot level [10,11].The task of processing this data while maintaining accuracy and meeting the time requirements becomes a real challenge.
Several methods have been proposed in the remote sensing literature to try to solve the automatic or semi-automatic parcel outlining problem.Most of them are based on image segmentation, edge detection algorithms, classification models, or combinations of these techniques.An object-based approach for extracting human-made objects, particularly agricultural fields from high-resolution images, was proposed in [12].This approach was able to extract regularly shaped objects by combining edge detection models with region-based segmentation.Da Costa et al. [13] developed an algorithm to outline vine plots automatically from very high resolution images by exploiting their textural properties.To differentiate between vine and non-vine pixels, they applied a thresholding method to the texture attributes of the image.In [14], a semi-automatic methodology for outlining field boundaries from satellite data was proposed.The authors first carried out segmentation using tonal and textural gradients, and the generated regions were then classified to obtain preliminary plot boundaries.Finally, they applied an active contour model [15] to refine the geometry of these boundaries.Turker et al. used perceptual grouping for automatically detecting sub-boundaries within existing agricultural fields from satellite imagery [16].This approach combined field boundaries and image data to carry out a field-based analysis.A Canny edge detector was used to detect the edge pixels, and then lines were identified using a graph-based vectorization method.
From the analysis of the state of the art, it can be concluded that approaches based on segmentation methods have several drawbacks including that (1) they are sensitive to intra-plot variability, which can result in the production of more segmentation than desired, and (2) most of these methods depend heavily on the correct selection of parameters (e.g., the similarity measured used to group image pixels), which needs prior knowledge of the landscape or tuning by trial and error.Moreover, variability in the sizes and shapes of the plots means that certain configuration parameters do not allow plots with the different characteristics needed for a landscape to be outlined properly.Approaches based on edge detection tend to produce more of the desired edges, mainly due to the presence of spatial patterns and image noise.The oversegmentation problem presented by the two approaches is directly related to the high spatial resolution of the images used in the segmentation process.Nevertheless, in the case of outlining plots, a very high spatial resolution (VHSR) is critical for managing very differently shaped and sized parcels in the same scene.Methodologies based on the superpixel (SP) concept have been proposed to deal with VHSR.These methodologies aim to reduce the influence of noise and intra-class spectral variability, preserving most edges of the images and improving the computational speed of later steps, such as classification, clustering, and segmentation [17].In fact, some approaches to outlining plots based on SPs are found in the literature.The authors in [18] combined a contour detection algorithm and simple linear iterative clustering (SLIC) to extract the cadastral boundaries from UAV orthoimages automatically.In [19], the problem of outlining plots was addressed as a machine learning problem.The authors first oversegmented a VHSR image, then labeled each pair of segments according to whether they belonged to an agricultural plot, and, finally, they trained a classifier using this information.The trained classifier was then used to segment the agricultural plots of other regions in the image automatically.The oversegmentation problem has been addressed in other areas, such as computer vision [20], video processing [21], and medicine [22], among others, by establishing a consensus between a set of different segments.This approach can also alleviate the parameter selection issue.Other works have addressed the problem of unsupervised parcel segmentation using time series.A procedure to identify different crops by combining information provided by an LPIS and a low spatial resolution image time series is presented in [23].However, LPISs are not always available; in fact, as has already mentioned, a suitable definition and update of an LPIS involve at least the semi-automatic outlining of the parcels.Nevertheless, some works in the literature that have applied time series for improving classification in agriculture.Thus, in [24], an operational crop identification strategy based on the use of multispectral and multitemporal signatures was proposed for classification at the parcel level.It proposes a combination of synthetic-aperture radar (SAR) and optical data registered on particular dates with the objective of satisfying crop temporal constraints.The authors proved that the use of SAR time series reduces the crop classification delivery time when the optical image is replaced by several SAR images.The integration of spectral and temporal features was proposed in [25] for annual cropland mapping.Even though the results presented in this paper showed that the methodology proposed is independent of in situ data and that it is capable of differentiating the croplands effectively, this methodology requires spatial baseline land cover information provided by different sources.A different crop classification approach is proposed in [26].Although an ensemble of multilayer perceptrons (MLPs) provides classified pixel-based and parcel-based maps from multitemporal satellite optical imagery, the labels of the training patterns were obtained through a ground survey.It should be noted that all these approaches are supervised, need additional information, and are not always available.To our knowledge, there is a lack of semi-or unsupervised methods that exploit temporal data for the purpose of outlining a high variety of parcels differing in size and shape.The use of temporal information to delineate agricultural plots appears promising since, during a growing season, plot boundaries in agricultural landscapes are relatively stable, while the phenological pattern of crops changes frequently [27].Some successful examples of combining superpixels and temporal information for change detection can be found in the literature [28,29].However, in this case, the parcel outlining problem could be seen as a special case of change detection, in which the pattern (parcel edges) tends to remain constant over time.
In this work, a new methodology based on consensus segmentation for outlining agricultural plots is presented.The proposed methodology combines segmentation on different scales-carried out using an SP method-with images registered on different dates to obtain a single segmentation of the agricultural plots.
This paper is organized as follows.Section 2 describes the imagery used and the ground truth built for the evaluation of the methodology proposed, which is explained in detail in Section 2.2.The results obtained are included and discussed in Section 3. The conclusions derived from this work are presented in Section 4.

Dataset
The study area comprises approximately 160 km 2 (16,000 ha) of fragmented agricultural plots located in the Lolol Valley, O'Higgins Region, Chile (Figure 1).The region is characterized by a temperate climate, with an agricultural season between September and April (the spring-summer season).The rainfalls are concentrated in the winter months (June-August).The landscape is characterized by very diverse sizes, ranging from small-scale farms on small plots-smallholdings-(<5 ha on average) to large legally constituted entities with medium and large plots (>50 ha).Three multispectral Plèiades-1 satellite images are available for analyzing the study area.The corresponding dates are 4 December 2017, 30 January 2018, and 25 February 2018.Table 1 summarizes the spectral and spatial characteristics of this imagery.The multispectral (MS) bands were geometrically corrected as well as co-registered.Moreover, the histograms of the images were adjusted to enhance the contrast.One scene containing mostly agricultural plots was selected and clipped from the three Plèiades images (as can be seen in Figure 2).The images under study correspond to a single agricultural season (2017-2018).As can be seen in Figure 2a-c, during an agricultural season, the main changes at the intra-plot scale correspond to the different phenological states of the crops.Furthermore, one ground truth map was obtained by manually outlining the agricultural plots in the clipped image for the date 30 January 2018.This map is used as ground truth to evaluate the result of the methodology proposed.Figure 2d sets out the polygons corresponding to the ground truth overlaid onto the color-composition image for the date 30 January 2018.

Methodology
The methodology proposed in this paper integrates segmentation at different scales, carried out by an SP approach, and different dates to obtain a final single segmentation of the agricultural plots.An overview of the proposed methodology is shown in Figure 3.
The first step of this methodology is the image segmentation at different scales.Segmentation was carried out by means of the multispectral SLIC algorithm proposed in [30].This method extends the original SLIC proposed in [31] by considering not only an RGB color space but a multispectral space.In this case, the clustering distance (Equation ( 1)) between two pixels p i and p j is composed of two components, one spatial (d s ) and the other spectral (d c ).
where x and y denote the position of the pixels.B represents the total number of bands in the multispectral image.The constants g and c influence the size of the superpixel and its compactness, respectively.The higher the g value, the bigger the superpixel; on the other hand, the bigger the c value, the more compact the superpixel.Segments at different scales are generated by fixing the value of the compactness parameter but varying the size of the SPs.Similar to [30], the number of SPs was selected to follow a dyadic progression.In this work, the three original images were segmented into 10 different scales using the SLIC algorithm.The parameters used in the segmentation are a compactness factor (c), which is 0.04 times the maximum spectral value contained in each of the image, and the number of SPs, which ranges from 2 8 (scale 1 ) to 2 17 (scale 10 ), on average.
The integration of the l segmentation for each date was carried out using a consensus process, which consists basically of a voting scheme that determines which pixels belong to the same region and which are part of the edges of objects in the image.
Thus, for each pair of adjacent pixels i and j, an E ij value is obtained by Equations ( 4) and ( 5), which takes into account whether the pixels belong to the same region.If the two adjacent pixels (the ith and jth) do not belong to the same region, that means they are part of the edges between these regions.E ij is defined as: where r k a and r k b represent the regions a and b, respectively, belonging to the kth segmentation, and S(r k a , r k b ) is an index of the similarity between the regions r k a and r k b .The larger the value of E ij , the stronger the separation between these regions.
The similarity index, analogous to that proposed by [32], is defined as a combination between the similarity indices of color and texture: where S color (r k a , r k b ) measures the color similarity.For each region r k u , a color histogram is obtained using 25 bins for each spectral band.Then, from the color histograms, a feature vector C u = {c 1 u , . . ., c H c u } with a length of H c = 25 × B is generated for region r k u (u can be a and b), where B is the number of bands of the multispectral images.The feature vector is normalized using the norm L 1 (also known as the Manhattan norm).Finally, the color similarity between the two regions r k a and r k b is calculated using the χ 2 statistic [33]: where S texture (r k a , r k b ) measures the texture similarity.Texture is calculated by means of steerable filters [34] using Gaussian derivatives in eight directions and σ = 1 as a basis.In this case, for each region r k u , a texture histogram is obtained using 10 bins for each spectral band and each filter direction, resulting in 80 features.Then, from the texture histograms, a feature vector T u = {t 1 u , . . ., t H t u } with a length H t = 80 × B is generated for each region r k u , where B is the number of bands in the image.The similarity in texture between two regions r k a and r k b is calculated as follows: The result of the consensus process is an edge map for each date that is generated with the information of the voting scheme, which is normalized to a range of 0-1.These maps combine different information, such as (i) different scales of segmentation through the different sizes of SPs, (ii) the dissimilarity in both texture and color between neighboring regions, and (iii) the probability of a pixel belonging to an edge.The closer to 1 the value of a pixel, the greater the probability of it belonging to an edge in the multiple segmentation scales, and, therefore, the contrast between regions separated by such a pixel at different scales is greater.
The next step in the methodology is the integration of the edge maps for all dates into the labeled boundaries (I edge ).In this work, this integration was carried out by averaging the n edge maps.
Low-pass filtering should be applied to reduce the noise present in I edge .In this work, a 3 × 3 median filter was used.
Since this averaged edge map generally includes open polygons, an ultrametric contour map (UCM) was calculated by means of the method proposed by [35].A UCM is an edge map with the remarkable property, i.e., it produces a set of closed curves when any threshold is set [36].The larger the threshold, the greater the contrast of the edges of the segments generated.The output of the UCM produces the outline of the final plots.In the end, a UCM is a soft representation of a segmentation that takes into account information from the edges of the image.The UCM has two inputs: the labeled boundaries (Equation ( 9)) and an image with boundary weights (IBW), which is defined by the Equation (10): (10) where SP scale 1 t 0 and SP scale 1 t n represent the SPs at scale 1 for the date t 0 and t n , respectively; and ∨ is the logical operator OR.A detailed description of the UCM algorithm can be found in [35].

Results and Discussion
In Figure 4, the SP segmentation obtained for the image registered on 30 January 2018 and corresponding to scale 4 and scale 8 is overlaid onto the color-composition image.It is possible to see the good adherence of the SPs to image objects in the three displayed cases.However, it can also be observed that the smaller the SP size (Figure 4a), the better its adherence to the edges, as well as the homogeneity of the pixels that compose it.While, in the segmentation scales with larger SPs (Figure 4a,b), the segments are less spectrally homogeneous, they still respect the borders between regions with a high contrast (the thickest edges).From the 10 segmentation scales for each date generated and consensus, the edge maps shown in Figure 5a-c were obtained.
As can be seen in Figure 5a-c, the edges of the plots tend to be more intense in all cases.In accordance with the methodology described, the next step is the integration of the three edge maps, one for each date, by calculating the average value.The average value is chosen as the basic form of integration because it allows for representing the information (borders) contained in all dates, thus reducing the effect of the appearance of borders in only one of the images, which is usually unusual for images of the same growing season (as mentioned in Section 1).In addition, the edges that appear as a result of phenological changes tend to decrease.The edge map obtained (Figure 5d) provides useful information for identifying the edges belonging to the plots; however, as mentioned above, it does not guarantee complete plots by applying a threshold, mainly due to edge-pixels not having the same value.Then, in order to obtain the final outlined agricultural plot map for the area considered, the UCM algorithm was applied to complete the edges.
Figure 6 shows the result of applying the UCM algorithm to the average of the three edge maps (Figure 5d).As can be observed, even though all edges have been closed, there is a lot of noise inside each plot.This noise is due to the fact that in the UCM, the totality of the probabilities of occurrence of the edges is represented.This is why it is necessary to determine a threshold value of the probability above which an edge is considered as such.Two indicators were used to determine this value objectively.The first was an error metric of the calculated edges with respect to a ground truth called boundary displacement error (BDE) [37,38]; the second metric considers the Shannon entropy, which represents the information content present in a border image: in this case, it is the number of edges present at the UCM.The BDE index measures the difference between two segmented images by averaging the displacement of boundary pixels.Specifically, the distance (d E ) between a boundary pixel (p s ) in the obtained boundary image (B s ) and the closest pixel (p gt ) in the ground-truth boundary image (B gt ) is used to define the error (disagreement) of each boundary pixel.The BDE index can be mathematically defined as in Equation (11).
where |B| represents the number of boundary pixels in image B, and d E is the Euclidean distance.The BDE index ranges within [0, ∞), where the lower its value, the better.The BDE index is plotted against the edge entropy in Figure 5d for each threshold value.As shown in Figure 7, the smallest BDE metric error occurs for threshold values that are close to 1, while the worst values occur for values close to 0. Conversely, entropy behavior presents low values of information content for high threshold values and high information content for high threshold values.For the purpose of determining a specific threshold value, a compromise was established between the BDE value and the entropy value.As can be observed, a threshold value of 0.5 provides the balance between the error and the amount of edge information.As can be seen, most of the edges agree between the two maps (Figure 8a,b), particularly in homogeneous areas with a high contrast between adjacent regions; however, discrepancies in plots with different tilled patterns can be also distinguished where the appearance of inner edges is present, as in plots within which anomalies are perceived due to poor agricultural practices, land heterogeneity, or crop diseases.As an example of the above, some areas where these changes occur are shown enclosed in circles.Consequently, only the edges stable over time-those that appear in the two images-remain in the final outlined plot map.

Conclusions
In this work, a new methodology based on consensus segmentation for outlining agricultural plots is presented.The methodology proposed combines segmentation at different scales-carried out using an SP method-and images registered at different dates to obtain a single segmentation of the agricultural plots.This methodology allows the outlining of agricultural plots of different sizes, shapes, colors, and textures.It is based on a consensus of segmentation (SP) at different scales and for different dates to determine the boundaries of the agricultural plots.The segmentation at different scales allows different-sized plots to be outlined.The use of SP for segmentation, which shows a good adherence to the edges for all scales, allows for differently shaped plots to be outlined.By highlighting the stable edges over time, the consensus of the segmentation reduces the intra-plot variability caused by the phenological stages present in a single growing season.The methodology also allows a threshold value to be determined in an objective way, which establishes a balance between the error and the amount of edge information.In particular, in this study, the threshold value determined for balance was 0.5.For lower threshold values, edges that are less stable over time appear on the outlined plot map, while for threshold values of up 0.5, only the more stable edges over time appear on the map.

Figure 2 .
Color compositions (3-2-1) of the Plèiades-1 imagery for three different dates and the ground truth of the study area.

E
ij (t) n < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 0 E f 9 M m X 0 H O D h g J V w r P r w E M 2 M R 4 = " > A A A C B n i c b V D J S g N B F H w T t x i 3 q E c v j U G I I G H G i 1 6 E g A g e I 5 g F k h h 6 O j 1 J a 8 9 C 9 x s h D A M e 9 T P 8 A T 2 J e v M D v H j 0 b+ w s B 0 2 s U 3 V V N b w q N 5 J C o 2 1 / W 5 m 5 + Y X F p e x y b m V 1 b X 0 j v 7 l V 0 2 G s G K + y U I a q 4 V L N p Q h 4 F Q V K 3 o g U p 7 4 r e d 2 9 O R 3 6 9 V u u t A i D S x x E v O 3 T X i A 8 w S g a q Z M v t H T s d x I 8 s d O r J E h J y 1 O U J W e d R F y n R d x P j Z Y z K b t k j 0 B m i T M h h b L z 9 X g H A J V O / r P V D V n s 8 w C Z p F o 3 H T v C d k I V C i Z 5 m m v F m k e U 3 d A e b x o a U J / r d j I q k 5 I 9 L 1 Q E + 5 y M 3 r + z C f W 1 H v i u y f g U + 3r a G 4 r / e c 0 Y v e N 2 I o I o R h 4 w E z G e F 0 u C I R l u Q r p C c Y Z y Y A h l S p g r C e t T M w S a 5 Y b 1 n e m y s 6 R 2 W H I M v z A 7 H M A Y W d i B X S i C A 0 d Q h n O o Q B U Y P M A z v M G 7 d W 8 9 W S / W 6 z i a s S Z / t u E P r I 8 f j 0 K b S A = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " G W 7 S 9 5 u o M 8 N M b F C A s C L G 9 S Q 4 D A 0 = " > A A A C B n i c b V A 7 S w N B G N y L r x h f U U u b x S B E k H B n o 4 0 Q E M E y g n l A c o a 9 z V 6 y Z u / B 7 n d C O K 5 W O / + C k F o r U T t / g I 2 l / 8 a 9 J I U m T j U 7 M w v f j B M K r s A 0 v 4 3 M 3 P z C 4 l J 2 O b e y u r a + k d / c q q k g k p R V a S A C 2 X C I Y o L 7 r A o c B G u E k h H P E a z u 9 E 9 T v 3 7 D p O K B f w m D k N k e 6 f r c 5 Z S A l t r 5 Q k t F X j u G E z O 5 i v 0 E t 1 x J a H z W j v l 1 U o T 9 R G s 5 n T J L 5 g h 4 l l g T U i h b X 4 + 3 w + F D p Z 3 / b H U C G n n M B y q I U k 3 L D M G O i Q R O B U t y r U i x k N A + 6 b K m p j 7 x m L L j U Z k E 7 7 m B x N B j e P T + n Y 2 J p 9 T A c 3 T G I 9 B T 0 1 4 q / u c 1 I 3 C P 7 Z j 7 Y Q T M p z q i P T c S G A K c b o I 7 X D I u 3 B l P x o v x O o 5 m j M m f b f Q H x s c P B E G d J w = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " G W 7 S 9 5 u o M 8 N M bF C A s C L G 9 S Q 4 D A 0 = " > A A A C B n i c b V A 7 S w N B G N y L r x h f U U u b x S B E k H B n o 4 0 Q E M E y g n l A c o a 9 z V 6 y Z u / B 7 n d C O K 5 W O / + C k F o r U T t / g I 2 l / 8 a 9 J I U m T j U 7 M w v f j B M K r s A 0 v 4 3 M 3 P z C 4 l J 2 O b e y u r a + k d / c q q k g k p R V a S A C 2 X C I Y o L 7 r A o c B G u E k h H P E a z u 9 E 9 T v 3 7 D p O K B f w m D k N k e 6 f r c 5 Z S A l t r 5 Q k t F X j u G E z O 5 i v 0 E t 1 x J a H z W j v l 1 U o T 9 R G s 5 n T J L 5 g h 4 l l g T U i h b X 4 + 3 w + F D p Z 3 / b H U C G n n M B y q I U k 3 L D M G O i Q R O B Ut y r U i x k N A + 6 b K m p j 7 x m L L j U Z k E 7 7 m B x N B j e P T + n Y 2 J p 9 T A c 3 T G I 9 B T 0 1 4 q / u c 1 I 3 C P 7 Z j 7 Y Q T M p z q i P T c S G A K c b o I 7 X D I

Figure 3 .
Figure 3. Overview of the methodology proposed.

Figure 4 .
Superpixels for different scales are overlaid onto the ground truth of the study area with red color.

(a) 4 5 .
December 2017.(b) 30 January 2018.(c) 25 February 2018.(d) Edge map, I edge Figure Edge maps for three dates and the edge map, I edge , obtained by averaging the three dates' edge maps.

Figure 6 .
Figure 6.Result of applying the ultrametric contour map (UCM) algorithm to the average of the three edge-maps.

Figure 7 .
Figure 7. Representation of error versus the edge entropy for different threshold values.

Figure 8
Figure 8 shows the final outlined plot maps obtained for two different threshold values overlaid onto the ground truth map.

Figure 8 .
Final outlined plot maps obtained for two different threshold values.

Table 1 .
Spectral and spatial characteristics of the multispectral image taken by the Plèiades-1 satellite.