Automatic Delineation of Forest Patches in Highly Fragmented Landscapes Using Coloured Point Clouds

Abstract

Accurate mapping of landscape features is key for natural resources management and planning. For this purpose, the use of high-resolution remote sensing data has become widespread and is increasingly freely available. However, mapping some target features, such as small forest patches, is still a challenge. Standard, easily replicable, and automatic methodologies to delineate such features are still missing. A common alternative to automated methods is manual delineation, but this is often too time and resource intensive. We developed a simple and automatic method from freely available aerial light detection and ranging (LiDAR) and aerial ortho-images that provide accurate land use mapping and overcome some of the aforementioned limitations. The input for the algorithm is a coloured point cloud, where multispectral information from the ortho-images is associated to each LiDAR point. From this, four-class segmentation and mapping were performed based on vegetation indices and the ground-elevation of the points. We tested the method in four areas in the north-western Iberian Peninsula and compared the results with existent cartography. The completeness and correctness of our algorithm ranging between 78% and 99% in most cases, and it allows for the delineation of very small patches that were previously underrepresented in the reference cartography.

1. Introduction

Accurate mapping and classification of landscape elements including different land covers (crops, pastures, forest areas, etc.) and ecosystems is key for land planning, natural resources management and biodiversity conservation [1,2,3]. In this sense, over the last few decades, a wide range of products and methodologies related to landscape mapping have been released (e.g., [2,4,5,6,7]). Available products for landscape delineation are mainly derived from remote sensors and cover a range of spatial, temporal and spectral resolutions that affect mapping accuracy [8,9]. For example, high temporal resolution data (i.e., one or more images per year) allow for monitoring ecosystems changes and dynamics (e.g., decline of forest ecosystems related with biotic or abiotic factors). High spectral resolution (i.e., availability of several spectral bands) allows the estimation of different types of indices (e.g., the widely used normalised difference vegetation index, NDVI) in order to support classification and mapping processes [10]. Finally, high spatial resolution images (i.e., pixel sizes smaller than 1–2 m [11]) allow for the delineation of small or very narrow features or patches. Small-sized forest patches, isolated trees, and non-forest woody formations are relevant elements of traditional landscapes of Europe [12]. Thus, their identification and delineation are frequently highlighted as key for the development of sustainable management and conservation tools for rural landscapes [13,14,15].

The manual delineation of features is the most obvious and direct method for mapping small-sized patches (e.g., [16,17]). For this reason, it was the basis for the first highly detailed landscape maps during the twentieth century. Manual delineation is often based on the use of aerial imagery (mainly ortho-images) with very high spatial resolution (typically with pixel sizes of 0.05–0.5 m [17]). However, these types of images are not always available in open repositories. The main limitation of these manual processes for patch delineation and identification is that they are very time consuming and subject to the interpretation of the operators.

Automatic and semi-automatic methods (i.e., typically unsupervised and supervised classifications, respectively) constitute an alternative to manual delineation and are commonly applied to satellite images from different sensors [10] such as Sentinel or Landsat. Although these data often provide higher levels of spectral and temporal resolution than the most common aerial images (i.e., aerial ortho-photographs [11,18]), they have, in general, lower levels of spatial detail. These limitations result in some small patches (or individual features) not being delineated.

Aerial light detection and ranging (LiDAR) provides a less time-consuming alternative to manual delineation. In addition, it provides higher spatial detail than traditional remote sensing imagery as well as information to characterise the spatial structure of the ecosystems (i.e., their height in relation to the ground floor level). Aerial light detection and ranging technologies are already well established and accepted as a useful tool for landscape mapping and classification processes (e.g., [7,19]). This type of dataset is becoming increasingly available. In many European countries, point clouds from their whole territory are available such as the National Air Ortho-photography Programme in Spain (PNOA) [20], Finland [21], Denmark [22] and Poland [23]. Other countries provide free point clouds of all their territories as is the case for the USA [24], Canada [25] and New Zealand [26]. In other world regions, including most of the South and Central Americas and Africa, data capture at the country level is undeveloped. However, several LiDAR studies [18,27,28] show the increasing interest of these countries in this data source. The data from these national LiDAR campaigns is usually provided with some sort of land cover/use classification. These classifications use point cloud geometry as a unique input and, in most cases, are based on a height normalization followed by multiple elevation thresholding and local texture analyses. References in the literature to these methods are numerous [29,30,31,32,33,34], and they often describe the use of existing software (e.g., Terrascan, LasTools, FME, Envi) or combinations of them. One of their main shortcomings is that they fail at differentiating covers with similar local textures and elevations.

In order to overcome some of the abovementioned limitations, our main aim was to develop a simple, automatic, and easily replicable method for mapping and classifying landscape patches using publicly available data: coloured LiDAR point clouds with multispectral data from aerial imagery. Four classes were mapped at very high spatial resolution (2 m pixel size): (i) forests, (ii) shrublands and low vegetation, (iii) buildings, and (iv) bare soil areas.

The developed method was tested in the highly heterogeneous landscape of the NW Iberian Peninsula (Galicia, Spain), and accurately identified different types of patches, particularly those which were woody dominated, included both forests and non-forest woody formations.

2. Materials and Methods

2.1. Study Area for Testing

To test our method, we selected four 2 × 2 km study areas located in Galicia, northwest of the Iberian Peninsula, European Atlantic Region (Figure 1). In general, the landscape of this area is highly heterogeneous with a large number of small forest patches combined with multiple types of land uses [35]. The predominant forest vegetation types are deciduous stands (mainly Quercus sp.) mixed with coniferous (Pinus sp.) and eucalypts (Eucalyptus sp.) plantations. On the one hand, the four areas are similar in terms of altitude (an altitude range of 50–250 m over sea level), type of terrain (slopes from moderate to high, between 5% and 35%), and ecological context (i.e., the same biogeographical regions). On the other hand, they show differences regarding landscape composition (e.g., higher values of forest cover in Area 4) and configuration (e.g., more fragmented patches in Area 1). See the results in Section 3.1 for further details.

2.2. Methodological Approach and Type of Data Used

Our algorithm used the LiDAR data available in the PNOA [20], which consist of LAS [36] binary files (the LAS format is used for managing LiDAR point cloud datasets) delivered in 2 × 2 km tiles of coloured point clouds. In addition to the information for the x, y, and z coordinates, every point in the cloud contains colour information from ortho-images in three bands: near infrared (NIr), red (R), and green (G). The data are freely available online in this format for the whole Spanish territory. The NIr sensor used in the data acquisition was a Leica RCD30 [37] which incorporates a monochrome charged coupled device (CCD) which is spectrally separated from the red, green, and blue (RGB) bands through a dichroic beam splitter device; thus, it can be considered as a near infrared-modified camera. In order to create distortion-free images of nominal focal length and pixel size, the NIr band was interpolated to the resolution of RGB. The point cloud and the spectral information were fused using a zenithal projection, so the colour attributes were added to each point in the cloud.

We developed a simple method that takes advantage of the two main aspects/attributes of the data we used: spectral and geometric information. As previously explained, the spectral data provided by coloured point clouds are modified for visualisation purposes and, therefore, not the most suitable for land cover classification. However, the robustness of some vegetation indexes (i.e., NDVI [38]) allows for a simple dichotomous classification (vegetated/non-vegetated land cover) even with this uncalibrated data. The calculation of the elevation of the elements/features on the ground from LiDAR point clouds has been extensively studied in the last two or three decades (e.g., [2,6]). Most of the existing methods generate a terrain model using regular grids or triangulated irregular networks. Once the terrain models are generated, calculating the elevation above them of every point in the original cloud is straightforward.

2.3. Methodology

2.3.1. NDVI-Based Segmentation

In the first stage of the method, NDVI is computed from NIr and R values [39] for each point. The procedure is implemented through a simple arithmetic operation, and it is performed on a matrix of size n × 2, n being the number of points in the cloud with two columns for R and NIr values. As a result, a new column/field is generated and added to the original point cloud that contains NDVI values, ranging from −1 to 1. Afterwards, points belonging to vegetated versus non-vegetated covers/features are separated by thresholding the NDVI values. This threshold may vary depending on the data source, study area, season, etc., but it is, in general, established in values around NDVI = 0 [40,41]. It is important to highlight that only calibrated NDVI values can be used to compare points in different locations, epochs or datasets; however, when the information is used to determine the distribution of relative variability (as it is in this study) the radiometric calibration of the data is not always indispensable. According to Xu et al. (2014) [41], several studies have obtained good results testing the potential of the near infrared-modified cameras (with ground-based, unmanned aerial vehicle UAV, or other airborne perspectives) to be used without calibration and transformation of the digital numbers (DNs) of the images into reflectance [11,42,43,44]. Figure 2 shows an RBG representation of the NIrRG values, and the result of the NDVI-based thresholding in a 2 × 2 km point cloud tile. For this study, we used a threshold = 0.0 for NDVI values.

2.3.2. Elevation-Based Segmentation

In the second stage of the method, the points from the cloud are height normalized, i.e., the height above the ground of each point is computed by modelling the terrain. The automatic generation of accurate terrain models from aerial point clouds has been tackled from many different perspectives [6]. In this study, we used the Cloth Simulation Algorithm (CSA), developed by Zhang et al. (2016) [45]. The CSA creates a regular grid that covers the whole extent of the point cloud. Then, the nodes of the grid are treated as a cloth with inverted gravity. The result is a cloth-like grid that adapts to the terrain, leaving out the elements that protrude from the terrain such as trees, buildings, or other volumetric structures. The resulting terrain grid (or virtual cloth) is smothered using a Laplacian filter. From this, distances of every point to the cloth are calculated (i.e., elevation of the points from the terrain model). The elevation is added as a new column/field to the original cloud. The same elevation threshold (3 m) was applied in the four test areas. This value represents the general minimum height of trees for this geographical area.

2.3.3. Segmentation and Classification

The original point cloud is subsequently segmented into four target classes. In order to perform this segmentation, threshold values for NDVI and elevation are considered in order to assign a class to each point (see Figure 2). The classification is performed as described in

Table 1. Descriptions of the four types of land uses used for our developed algorithm and for the two reference maps used for validation: Corine Land Cover 2018 (CLC18) [46] and Spanish Forest Map 1:25,000 (MFE25) [47]. In the third column, the numbers in brackets show the code from CLC18 classification.

Class	Algorithm	CLC18 (Corine Land Cover 2018) Classes Included	SFM25 (Spanish Forest Map 1:25,000) Classes Included
Forest and trees	High (i.e., above the elevation threshold) and vegetated (i.e., the normalised difference vegetation index NDVI threshold)	Forests (31)	Polygons where tree cover was >50% (or where tree cover > shrub/grass cover)
Shrub and low vegetation	Low (i.e., below the elevation threshold) and Vegetated	Artificial vegetated areas (14) Permanent crops (22) Pastures (23) Heterogeneous agricultural areas (24) Shrub and/or herbaceous vegetation (32)	Areas where shrub/grass cover was >50%
Buildings and artificial areas	High (i.e., above the elevation threshold) and non-vegetated (i.e., below the NDVI threshold)	Urban fabrics (11) Industrial areas (12)	No applicable (this type of cover was not included in this dataset)
Bare soil	Low (i.e., below the elevation threshold) and non-vegetated (i.e., below the NDVI threshold)	Arable lands (21) Mines and construction sites (13) Open spaces with little vegetation (33)	Areas without vegetation cover

: (i) points are tagged as “Forest” if they are high (above the elevation threshold) and vegetated (above the NDVI threshold); (ii) “Shrubland or low vegetation” if they are low (below the elevation threshold) and vegetated (above the NDVI threshold); (iii) “Building” if they are high (above the elevation threshold) and non-vegetated (below the NDVI threshold); and (iv) “Bare soil” if they are low (below the elevation threshold) and non-vegetated (below the NDVI threshold). The final result (classified point cloud) is transformed into a raster. The ground sample distance (GSD) of the raster may vary depending on the density of the point cloud, so there is at least one point at the majority (or all) of the positions of the grid/raster (pixels). In cases where several points are at the same position of the grid, a preference rule applies in the following order: (i) forest, (ii) buildings, (iii) shrubland and low vegetation, and (iv) bare soil. In this way, if there are points from different classes in the same pixel, the pixel is tagged as the point of highest order (e.g., if there are “forest” points, the pixel is tagged as so; if not, if there are “building” points, the pixel is tagged as “building”, and so on). If there is a pixel with no points in it, it will be tagged as the majority of its neighbours are in order to fill the possible gaps.

2.4. Validation

We validated our algorithm using a regular grid of 32 × 32 points (n = 1024; distance between consecutive points: 60 m; see Validation grid in Figure 2) for each study area. The class values for each point (ground truth) were extracted from visual inspection of the coloured point clouds by two different operators. We also compared the performances of the algorithm for the four 2 × 2 km study areas with reference cartography from two different reference maps:

(i)

Corine Land Cover (CLC18) [46]. This version was updated in 2018. It provides a map of 44 land cover hierarchical classes for all of the European Union with a minimum mapping unit (MMU) of 25 ha and a minimum width of linear elements of 100 m;

(ii)

Spanish Forest Map 1:25,000 (SFM25) [47]. This thematic map focuses on forest ecosystems, providing an accurate delineation of tree-dominated ecosystems with multiple applications for forest and land planning at different scales (e.g., [48,49]). The SFM25 has a higher spatial resolution than CLC18 with an optimum pixel size of 10–15 m (according to recommendations in Hengl (2006) [50]) and a MMU ranging from 0.5 to 2 ha). For each landscape polygon delineated, it includes the percentage of tree and shrub cover and the three dominant tree species where these exist.

To compare the three different datasets, we clustered different classes (both in CLC18 and in SFM25) to create new broader classes comparable to those generated by our algorithm. We created a typology of four main classes based on the existence of vegetation cover and the height of the landscape elements. Then, we adapted the legend of the two reference maps to these four classes. For example, some semi-natural areas from CLC18 (i.e., shrublands) were separated from forests and combined with other types of low vegetation. Further information is provided in Table 1. Similarly, we re-classified SFM25 using the information from each individual polygon into the four classes previously mentioned (Table 1).

The performance of the method was measured by calculating confusion matrices (a widely accepted method for performance estimation [51]) from the values in the validation grid and the three datasets: results from the algorithm, CLC18 and SFM25. We computed the completeness and correctness for each class (forest/trees, shrubland, buildings, and bare soil), area, and dataset. Finally, we estimated the total accuracy and the Cohen’s kappa coefficient for our algorithm and for the two reference maps.

Completeness (CM; Equation (1)) provided information about the detection rate, i.e., the proportion of points that were visually assigned to a specific class in the validation grid and were correctly detected.

$$nCorr=number of points in the validation grid correctly assigned to a class$$

$$nReal=total number of points in the validation grid belonging to a class$$

$$nPred=total number of points in the validation grid assigned to a class$$

$$CM=\frac{nCorr}{nReal}$$

(1)

Correctness (CR; Equation (2)) provided information about the correct classification rate, i.e., the proportion of points from the validation grid that were correctly assigned to the class they belong to.

$$CR=\frac{nCorr}{nPred}$$

(2)

Accuracy (AC; Equation (3)) provided information about the overall completeness of the method, i.e., the proportion of the total number of points that were correctly assigned to their class.

$$AC=\frac{{\displaystyle \sum }nCorr}{{\displaystyle \sum }nPred}$$

(3)

Cohen’s kappa coefficient (CK; Equation (4), according to References [52,53,54]) provided information about accuracy, giving compensation to possible random hits.

$$CK=\frac{{\displaystyle \sum }nPred \times {\displaystyle \sum }nCorr- {\displaystyle \sum }\left(nReal \times nPred\right)}{\left({\displaystyle \sum }nPred\right)²- {\displaystyle \sum }\left(nReal \times nPred\right)}$$

(4)

3. Results

3.1. Delineation of Landscape Elements

The method proposed here allowed the automatic mapping of the four study areas at high spatial resolution (pixel size 2 × 2 m). The classification of the study areas, which covered 4 km², only required seconds of data processing using a standard computer. In general, the four study areas show that low vegetation (including shrubs, pastures, and crops) were the dominant cover (40–70%; Figure 3). Forest covered the largest areas with different levels of fragmentation (e.g., larger and more continuous patches in Area 2 and higher fragmentation in Area 4) or shapes (predominant lineal patches in Area 1). Tree patches were scarcer in Area 1 (<20%), where lineal shapes were predominant which represented forested hedgerows and similar structures. Areas covered by buildings and other artificial elements were identified within the four areas but, in general, showed lower cover values.

The delineation derived from the two reference maps (i.e., SFM25 and CLC18) provided different results to those obtained by our algorithm. In general, both datasets showed lower levels of occurrence for areas covered by forest patches in our study cases (lower than 45% except for Area 2). Comparing our findings with SFM25, the delineation of forested areas was similar, although a higher number of linear elements covered by woody vegetation were identified by our algorithm (see for example Area 4 in Figure 3). Our methodology also showed a higher level of performance for detecting woody hedgerows in comparison with CLC18. In fact, the last did not show any forest/tree/vegetation cover for Area 1, meanwhile our algorithm showed a tree cover of <20%. In addition, and in contrast with our algorithm, the two reference maps showed that those areas with bare soil (mainly arable lands) constitute the dominant class (40–70%; Figure 3). Finally, artificial covers (e.g., roads or paved areas) and buildings were practically absent in our study areas (Figure 3), but our algorithm detected existing ones, whereas CLC18 did not detect any (except for some pixels from Area 3), and SFM25 did not include this type of use.

3.2. Validation of Our Results and Comparison with Other Datasets

The results obtained with our algorithm showed high levels of performance based on the four different metrics (i.e., completeness, correctness, accuracy, and Cohen’s kappa) previously described (Figure 4). Our algorithm showed high values of accuracy (>92.8%) and Cohen’s kappa (>87.2%) for the four study areas. These metrics showed a lower performance for both SFM25 (<77.4% for AC and <57.9% for CK) and CLC18 (<69.4% for AC and <41.0% for CK). In relation to the areas defined as forest patches, at least 97.9% (CR) of the points were correctly classified. In addition, at least 89.8% (CM) of the points that were forests were detected. Detection and classification of the forest class showed the highest levels of performance among the four types of uses included in our study. The other types of land use also showed high levels of accuracy when using our own algorithm, especially the low vegetation class with values higher than 95.0% for completeness and 88.4% correctness. On the contrary, the two reference maps showed lower levels of accuracy for both completeness and correctness metrics (Figure 4). For CLC18, completeness of forest cover ranged from 0.0% to 70.3%, and correctness from 0.0% to 79.0%. These metrics showed higher levels of performance of SFM25 for the forest patches class: completeness 26.5–85.2% and correctness: 70.4–87.7%. Regarding the other vegetated class (shrublands and low vegetation), our algorithm also showed a better performance (correctness ranging from 88.4% to 96.6% and completeness from 95.0% to 99.6%) than the two reference maps (CLC18: 52.1–69.4% for correctness and 73.3–100% for completeness; SFM25 63.3–72.6% for correctness and 76.3–97.3% for completeness). The remaining two classes showed lower levels of performance for our algorithm (bare soil CM: 77.8–100% and CL 84.8–90.0%; and buildings CM: 50.0–100% and CL: 89.1–98.8%), but they were rarely identified by the two reference maps. As a consequence, validation metrics for those two classes were not calculated for CLC18 and SMF25.

4. Discussion

The algorithm developed in this study allows for the identification and delineation of different types of landscape features and land uses at the local level including forest/tree/woody vegetation patches of small size (<0.5 ha). This automatic method, based on the combination of point clouds derived from LiDAR data and aerial imagery, produces a higher spatial resolution output in comparison with other available datasets such as satellite imagery. The information used is freely available for the whole area of Spain, and every point within the point clouds contains spectral information in three bands (G, R and NIr) with different potential applications for classification of landscape patches by using vegetation indexes [10].

Spain is not the only country which has freely available LiDAR data of its whole territory. In recent years, a rising number of countries have conducted data collection for this kind of product. The increasing range of applications in different fields and the reduction of acquisition costs appear to be the main factors which have contributed to this rise. Moreover, some countries are investing efforts not only in the acquisition of data for their whole area but also in its periodic update [18,20,21,22,23,24,25,26,27,28]. Therefore, there is an increasing body of LiDAR datasets which will allow the use of the developed algorithm in many areas worldwide. In countries where the data are updated periodically, the algorithm can also be used for analysis of landscape dynamics and changes.

The spectral data from the aerial photographs used for the analysis were acquired by a non-calibrated digital camera that captures near infrared and visible wavelengths. These kinds of consumer-grade digital cameras are now established as a cost-effective method of monitoring plant health and phenology [38]. The combined use of infrared and visible observations provides more information than the use of just visible bands due to the strong contrast between infrared and visible reflection of vegetation, the high correlation of the three visible bands, and the possibilities of using spectral indices such as the NDVI. However, rigorous radiometric corrections are not applicable for images that have been processed and optimised for image visualization. In these cases, the pixel value may not represent the originally measured reflected brightness values and, therefore, it would not be possible to calculate the reflectance.

Our algorithm provided good results for identifying and delineating forest and woody patches (performance evaluated by confusion matrices and correctness and completeness values). In addition to this, our algorithm identified a higher level of vegetation cover than the other two alternative datasets compared (i.e., CLC18 and SFM25). Some of the differences, most of which related to the “low vegetation areas”, may be derived from areas of crops with temporal cover (arable areas) and showed positive NDVI values. This was when the aerial images were taken when the areas were covered by annual crops. Our methodology, despite having different and relevant strengths, is limited by the availability of adequate LiDAR and imagery information. Another shortcoming of our method is that it cannot establish the differences among forest types in terms of species composition, because this type of analysis demands the use of additional spectral data [55]. In addition, regarding the aboveground elevation, our method involves the generation of accurate terrain models from aerial point clouds. The generation of this type of model is not easy and has previously been tackled from many different perspectives [6]. The selection of the best terrain modelling algorithm remains subjective, however, as no single algorithm is objectively more accurate in all possible scenarios [56].

Finally, there are additional available sources of data that could be useful for similar aims to those assessed in this study. The European Union’s Earth Observation Programme “Copernicus” provides meaningful data. Copernicus supplies information with intermediate spatial resolution (pixel sizes between 20 and 100 m) which is periodically updated and includes relevant descriptors such as tree cover density or dominant leaf type [57]. Comparison between our algorithm and Copernicus-derived information was not included in this study, because Copernicus does not provide information about other woody elements (i.e., shrubs), bare soil areas, etc. In addition, the spatial resolution of the products derived from Copernicus can be limiting for delineating small forest patches, such as lineal woody elements, that are often narrower than 10 m [35]. In addition, additional applications can be derived from our methodology. For example, the digital model with the height above ground of the different vegetation covers can be used for multiple aims related to forest and land planning (e.g., estimation of biomass stocks [58] or fuel modelling for studies related to fire behaviour [59,60].

5. Conclusions

This study presented an automatic and easily reproducible algorithm for the classification and mapping of landscape patches with high spatial resolution (an approximately 2 m pixel size). The methodology developed combined information from LiDAR digital elevation models (for calculating the height of the different vegetation and artificial elements) and spectral information for G, R, and NIr bands derived from aerial imagery (for identifying the occurrence of vegetation covers). An important aspect of the proposed methodology was that we only used information available in public repositories.

The algorithm was tested on four different highly heterogeneous and highly fragmented areas in NW Spain, showing higher levels of accuracy than the two reference maps used as comparisons. The methodology presented here can constitute the basis for future high-detailed assessments at the landscape level of ecosystem dynamics and land cover changes. The method is easily reproducible, and for a 4 km² area, the application only required seconds of data processing on a standard computer. The outputs generated were small-sized files (approximately 4 MB) which can allow the characterisation of larger study areas.