Forest Cover Mapping Based on a Combination of Aerial Images and Sentinel-2 Satellite Data Compared to National Forest Inventory Data

Abstract

Research Highlights: This study developed the first remote sensing-based forest cover map of Baden-Württemberg, Germany, in a very high level of detail. Background and Objectives: As available global or pan-European forest maps have a low level of detail and the forest definition is not considered, administrative data are often oversimplified or out of date. Consequently, there is an important need for spatio-temporally explicit forest maps. The main objective of the present study was to generate a forest cover map of Baden-Württemberg, taking the German forest definition into account. Furthermore, we compared the results to NFI data; incongruences were categorized and quantified. Materials and Methods: We used a multisensory approach involving both aerial images and Sentinel-2 data. The applied methods are almost completely automated and therefore suitable for area-wide forest mapping. Results: According to our results, approximately 37.12% of the state is covered by forest, which agrees very well with the results of the NFI report (37.26% ± 0.44%). We showed that the forest cover map could be derived by aerial images and Sentinel-2 data including various data acquisition conditions and settings. Comparisons between the forest cover map and 34,429 NFI plots resulted in a spatial agreement of 95.21% overall. We identified four reasons for incongruences: (a) edge effects at forest borders (2.08%), (b) different forest definitions since NFI does not specify minimum tree height (2.04%), (c) land cover does not match land use (0.66%) and (d) errors in the forest cover layer (0.01%). Conclusions: The introduced approach is a valuable technique for mapping forest cover in a high level of detail. The developed forest cover map is frequently updated and thus can be used for monitoring purposes and for assisting a wide range of forest science, biodiversity or climate change-related studies.

1. Introduction

Forest cover is a key variable of interest to science, management and reporting [1]. Up-to-date and accurate information on the dynamics of forest cover is essential for the conservation and management of forests [2,3,4], carbon accounting efforts and the parameterization of a broad range of biogeochemical, hydrological, biodiversity and climate models [5]. A forest cover map provides information on the size, shape and spatial distribution of forests and thus assists in classifying the landscape into patterns. The spatial distribution of these landscape patterns is linked to the ecological functionality of the landscape [6] and provides new perspectives for ecological connectivity studies [7]. Therefore, the assessment of forest cover is aimed at facilitating decisions on biodiversity conservation and reforestation programs [8]. Furthermore, forest maps are crucial for global environmental change assessment and local forest management planning [7]. Especially in a changing climate, forest structures are changing very rapidly, meaning that knowledge of stocked areas and their changes is very important. Accordingly, forest information must be accurate, spatially detailed and up to date [9]. Due to the aforementioned reasons, there is an important need for a forest cover map containing a high level of detail.

Forests are inventoried for operational planning and forest management [9]. In many countries, National Forest Inventories (NFIs) are based on field surveys [10] and collect information on national forest resources. This enables strategic planning and policy development at the national level [9]. However, the aforementioned demand for information about forest cover exceeds the scope and design of most of the existing forest inventories [11]. Traditionally, forest cover maps are produced by visual interpretation of aerial images in combination with field surveys. Consequently, their development is time-consuming and therefore limited to relatively small areas [7].

In the last decade, various remote sensing-based methodologies have been used to map forest cover. The scales range from regional [12,13,14,15,16,17,18,19,20,21,22,23], nation-wide ([7] for Switzerland, [1] for Canada, [24] for Kyrgyzstan, [25] for Norway, [26] for Canada’s Boreal Zone, [27] for Ukraine plus Eastern Europe and [3] for Japan), pan-European [28,29,30] and post-Soviet Central Asian [8] to the global level [31,32,33,34,35,36]. Besides aerial images and light detection and ranging (LiDAR) data, mainly satellite data have been used for the mapping and monitoring of forested areas and vegetation. Most studies involved either Sentinel-2 data [12,13,17,18,20,25,29] or LANDSAT data [1,26,27,30,31,36,37]. In contrast, the authors of [15] used SPOT data for forest cover and forest degradation mapping, whereas other studies used radar data like TanDEM-X [33], airborne S-Band radar [16] or Sentinel-1 [14]. The studies of [7,23], both based on digital surface models (DSMs) from image-based point clouds, produced wall-to-wall forest maps [7] and modeled fractional shrub/tree cover [23]. Furthermore, the work of [35] characterized global forest canopy cover using spaceborne LiDAR. In this context, the authors of [19] classified forest and non-forest based on full waveform LiDAR. Besides single sensor approaches using only one type of satellite or airborne data, multisensory approaches have been used for forest cover mapping. For example, [38] used LiDAR in combination with aerial images, [28] used SPOT 4/5, Corine Land Cover 2006 and LISS III, [3] used LANDSAT and Sentinel-1, [8] used LANDSAT together with MODIS data and [24] used a hybrid fusion strategy from Globe Land 30 2010 and USGS tree cover 2010. The authors of the aforementioned studies were mainly driven by the motivation to conduct change studies. Using time series analysis of LANDSAT images, the authors of [31] characterized global forest extent and change from 2000 to 2019. The studies of [12,13,20,21] each focused on the analysis of forest succession detection, height differences and forest cover dynamics as well as on forest conservation and deforestation. For validation, stand-level observations [25], maps of forest structure [26], high-resolution satellite data [15], LiDAR data [14], Bing or Google maps [3,37], Google Earth [18], orthoimages and terrestrial surveys [19] or NFI data [1,7,18] have been used.

Although there are numerous approaches in the literature to derive forest cover and forest cover changes, the derivation of forest areas and thus the differentiation between forest and non-forest is still a challenge. While available global or pan-European forest maps have a low level of detail and forest definition is not queried, local administrative data are, in many cases, oversimplified or out of date [7]. Another challenge is posed by the different definitions of the term “forest” at the regional, national and international levels. The term “forest” is a summarizing generic term which includes various forms of land cover, including tree-free forms, such as open spaces, wild meadows or traffic routes. In contrast, there are also tree-covered landscapes which are not defined or perceived as forests, including field shrubs, Christmas tree plantations or short-rotation plantations [39]. Land use is a key parameter in forest definitions but hardly assessable with remotely sensed data [7]. The most important quantitative characteristics of forest areas include minimum size, minimum width, tree height, minimum stocking level and stocking duration. A derivation of forest from remote sensing data according to the abovementioned definitions is not directly possible. However, it is technically feasible to record stocked and unstocked areas and to determine measurable criteria like height, area size and stand width [39].

The main objective of the present study was to generate a spatio-temporally explicit forest cover map, which takes the German forest definition into account. The key criteria are (1) minimum tree height, (2) minimum crown closure, (3) minimum width and (4) minimum area. The forest cover map based on the results of this study comprises the whole region of Baden-Württemberg, Germany, covering around 35,751.5 km². Furthermore, the forest cover map was compared to NFI data with respect to the presence or absence of forest. Finally, the differences between the forest cover map and NFI data were analyzed, categorized and quantified. To generate the forest cover map, we used a multisensory approach involving both aerial images and Sentinel-2 data. In fact, Sentinel-2 is designed for vegetation monitoring [40] and therefore multiple studies and projects exploited these data to generate vegetation maps. The spatial, temporal and radiometric resolution of Sentinel-2 enabled tree cover mapping with a medium level of detail; however, in combination with orthoimages from aerial image flights, a high level of detail could be achieved. With digital elevation models from image-based point clouds, vegetation height [39] was queried, resulting in the production of a tree cover and forest cover map. The introduced forest cover map does not include information about land use. In this context, biodiversity-relevant conditions or changes such as storm damages, permanent or temporary forest gaps or forested areas outside the forest area are mapped.

2. Materials and Methods

2.1. Study Area

The study area comprises the federal state Baden-Württemberg, located in southwest Germany. Baden-Württemberg is located between 7°30′ and 10°30′ E and 47°32′ and 49°47′ N. It covers 35,751.5 km², with an altitudinal range between 84 and 1493 m a.s.l. The state is characterized by different types of landscapes and is divided into seven main growing regions (see Figure 1).

Baden-Württemberg has a varying regional proportion of forest cover with an average coverage of approximately 37.26% (±0.44%), which equals 13,319.6 km² (±157.8 km²) in total. The forest ownership types are distributed as follows: corporate forests account for 40% of the forest area, private forests come to 35.9% and state-owned forests have 23.6%, whereas federal forests only cover 0.5%. Five tree species, namely, Norway spruce (Picea abies (L.) H. Karst), beech (Fagus sylvatica L.), silver fir (Abies alba Mill.), pine (Pinus sylvestris L.) and oak (Quercus sp. L.), occupy over 75% of the forest area and determine the forest character. Smaller percentages are found for ash (Fraxinus excelsior L.), Sycamore maple (Acer pseudoplatanus L.), Douglas fir (Pseudotsuga menziesii (Mirb.) Franco), larch (Larix decidua Mill.), hornbeam (Carpinus betulus L.), birch (Betula sp. L.) and alder (Alnus sp. Mill.). A total of 50 different tree species have been recorded in Baden-Württemberg, most of them with minimal proportions [41]. More detailed information on forests of Baden-Württemberg is provided in the third NFI report [41].

2.2. Remote Sensing Data

This section describes the used remotely sensed data. Aerial orthoimages and normalized digital surface models (nDSMs) allowed the extraction of vegetation and vegetation heights in a high level of detail, whereas Sentinel-2 data were used to derive vegetation in an additional, generalized approach.

2.2.1. Aerial Images

The orthoimages and nDSMs used in this study were processed from aerial images acquired by airborne image flights between 2011 and 2014. The aerial images provided by the state agency of spatial information and rural development of Baden-Württemberg (LGL) [42] contained the four bands red, blue, green and infrared (RGBI) at a radiometric resolution of 16 bits. Flight conditions such as time of day, season and weather conditions, as well as flight settings such as front/side overlap and camera type, as well as image-matching parameters such as the version of the image-matching software (SURE) and ground sampling distance (GSD), varied between the flight missions. The flights were conducted in the vegetation period, with the earliest flight on 15 April and the latest flight on 19 October. Front overlap was between 59% and 70%, whereas side overlap was between 28% and 40%. The aerial images were taken with UltraCam Eagle, UltraCam Xp, DMC 01 or DMC II. Photogrammetric point clouds and orthoimages were processed with a GSD of 40 or 50 cm using the software SURE of nFrames [43]. Detailed information on the settings of the image flights can be found in

Table A1. Settings of the image flights.

Year	Mission	Date of Acquisition	Overlap [%]	Camera	SURE Version	GSD [m]
2013	1	12 July 2013	60/33	UltraCam Eagle	1.3.1.132	0.4
	2	13 July 2013 14 July 2013	62/37	UltraCam Eagle	1.3.1.336	0.5
	3	17 June 2013	64/38	UltraCam Xp	1.3.1.336	0.5
	4	19 June 2013	59/28	UltraCam Xp	2.2.1.1189	0.5
	5	17 June 2013	63/46	UltraCam Eagle	1.2.1.209	0.4
	6	18 June 2013 14 July 2013	62/40	UltraCam Xp	1.2.1.209	0.4
	7	13 July 2013 16 July 2013	64/38	UltraCam Xp	2.2.1.1189	0.5
	8	6 June 2013 7 June 2013 8 June 2013	63/34	UltraCam Eagle	2.2.1.1189	0.5
	9	19 June 2013 12 July 2013	59/31	UltraCam Xp	2.2.1.1189	0.5
2012	1	24 June 2012 30 June 2012	59/29	UltraCam Xp	3.0.0.0	0.5
	2	25 May 2012 23 July 2012	59/32	UltraCam Eagle	2.3.1.66	0.5
	3	23 July 2012	60/30	UltraCam Eagle	3.0.0.0	0.5
	4	15 April 2012	59/33	UltraCam Xp	2.3.1.66	0.5
	5	1 August 2012	61/40	DMC II	2.2.1.1212	0.5
	6	26 May 2012 25 July 2012 27 July 2012 1 August 2012 19 October 2012	61/39	DMC II	2.3.1.66	0.5
	7	23 July 2012	59/34	DMC 01	2.3.0.38	0.5
	8	25 July 2012	61/43	UltraCam Xp	1.3.1.336	0.4
2011	1	5 May 2011 6 May 2011	66/34	DMC 01	3.0.0.0	0.5
	2	6 May 2011	62/32	UltraCam Xp	3.0.0.0	0.5
	3	6 September 2011 10 September 2011 16 September 2011	61/32	UltraCam Xp	3.0.0.0	0.5
	4	3 September 2011	62/31	UltraCam Xp	3.0.0.0	0.5
	5	5 May 2011	64/31	UltraCam Xp	2.3.0.38	0.5
2014	5	26 September 2014	61/34	UltraCam Xp	1.3.1.284	0.4
	6	22 June 2014 26 June 2014	70/45	UltraCam Eagle	1.3.1.284	0.5
	7	6 June 2014	60/31	UltraCam Xp	1.3.1.392	0.5

. As shown in Figure 2, the flight missions covered 99.7% of the study site.

2.2.2. Digital Surface Models

DSMs generated from photogrammetric point clouds served as the basis for deriving forest heights. Therefore, a combination of lasgrid and las2tin from LASTools [44] was applied. With lasgrid, the Z-value of the highest point per pixel was returned. Voids were filled with a square search radius of 3 pixels. Remaining voids were interpolated using the las2tin algorithm to triangulate the point cloud. The geometric resolution was set to 1 m. In order to obtain forest heights, the difference between the DSM value and the corresponding terrain height (=nDSM) was calculated. For this purpose, we used the LiDAR-based digital terrain model (DTM) of the LGL with a GSD of 1 m [45]. A detailed description about the applied methods for deriving aerial image-based nDSMs can be found in [46].

2.2.3. Sentinel-2 Data

In the frame of the Copernicus program, the European Space Agency (ESA) has launched the Sentinel-2 optical imaging mission. The Sentinel-2 mission involves a constellation of two polar orbiting satellites: Sentinel-2A and Sentinel-2B. Each satellite is equipped with the optical imaging sensor Multi-Spectral Instrument (MSI). The Sentinel-2 satellites provide images in the visible and infrared spectrum between 443 and 2190 nm and therefore are optimized for land surface observation [47]. Sentinel-2A was launched in 2015 and Sentinel-2B followed in 2017 [48]. According to [49], the resolution of up to 10 m and the scanning width of 290 km are ideal for detecting changes in vegetation and for making harvest forecasts, mapping forest stands or determining the growth of wild and agricultural plants.

In this study, we used Level-2A products of Sentinel-2. The processing of Level-2A products was conducted with Sen2Cor [48] and consisted mainly of scene classification and atmospheric correction. Scene classification is a necessary input for atmospheric correction. Atmospheric correction is performed in order to obtain bottom of atmosphere (BOA)-corrected transforms of the multispectral satellite imagery [50]. More information can be found in [48,50]. The image collections of Sentinel-2 (Level-2A) were compiled for the growing seasons (1 May to 31 September) of the years 2017, 2018 and 2019 taking all images into account with a cloud cover of < 25%. The time frame (2017 to 2019) was chosen due to the availability of Sentinel-2 (Level-2A) data, which have been available since 2017. The time gap between both NFI data and aerial images to the Sentinel-2 data was not expected to have much influence on the resulting forest cover map, as the vegetation layers were mainly used to remove buildings. The forest cover map is mainly determined by vegetation height derived via nDSM. To cover Baden-Württemberg, 9 Sentinel-2 granules were required. The used granules and the number of scenes per year are shown in

Table 1. Number of Sentinel-2 granules used to calculate a vegetation mask.

Sentinel-2 Granules	2017	2018	2019	∑
32UMU	3	4	6	13
32UMV	4	5	7	16
32UNA	8	8	9	25
32UNV	10	5	11	26
32UNU	5	1	7	13
32TMT	9	5	15	29
32TLT	14	19	36	69
32TNT	6	1	6	13
32ULU	7	19	34	60

. Since the date of recording and cloud cover varied between the Sentinel-2 granules, the number of the used scenes per granule and year differed.

2.3. NFI Data

The German NFI records a forest in its size and structure and thus provides information on the state of the forest and its development. In Germany, the first NFI was carried out in the years 1986 to 1988. The second NFI followed from 2001 to 2002 and the third NFI was from 2011 to 2012. Since 2010, a ten-year cycle has been in force for the NFI [41].

The present study used plot data of the latest NFI (2011–2012). The sampling design of the inventory utilizes a 2 × 2 km grid aligned on a systematic, national sampling grid. Grid points define the lower left corners of a square inventory tract of 150 × 150 m, whereas corner points of the tracts mark the centers of the NFI sampling plots [51]. Plots located in the forest were geolocalized with the Global Navigation Satellite System (GNSS) device MXbox of GEOSat GmbH [52] including a correction signal. At each tract corner point located in the forest, nested circular sub-plots of 5 different radii and a Bitterlich sub-plot [53], so-called angle count sampling (ACS), were established to record the set of inventory variables regarding trees, stand structure and site characteristics. With ACS, sample trees are selected with a probability proportional to their basal area, meaning that the range of each plot differs [54].

With the help of administrative forest cover maps and aerial photographs, new, previously non-forest tract corners in forests are identified (“preclarification”). Whether a tract corner is declared as being inside or outside of the forest is finally decided in the field (“forest decision”, see Figure 3). The coordinates of the visited forest plots were corrected after the best attempt at reaching the point, whereas the non-forest plots have a tentative plot location. Out of the latest NFI, 35,744 NFI plots were available, while 13,299 were covered by forest, 425 were temporary or permanent unstocked areas and 22,020 were outside of the forest areas. Amongst other variables, tree heights were measured for a subset of sample trees. Subsequently, the tree heights of the remaining sample trees were modeled. Furthermore, information about land use, forest edges and accessibility was recorded.

To analyze forest cover within an NFI plot, the NFI tract corners were buffered with a radius of 10 m. This corresponds to the median distance between the NFI plot center and the maximum measured distance of all corresponding sample trees. The authors of [54] considered this approach as reasonable. In a further step, the percentage of forest cover within each plot was calculated. Since forest cover layer refers to land cover and NFI data to land use, the NFI data were filtered based on the available NFI plot information. For this purpose, only those NFI plots which were located in a forest with trees higher than 5 m or outside the forest area were evaluated. As a consequence, the following NFI plots were removed subsequently:

• NFI plots at non-wood areas (n = 425);
• NFI plots where the forest cover map could not be calculated due to clouds or artefacts (n = 153);
• NFI plots which were not accessible (n = 78) and therefore no forest inventory was conducted;
• NFI plots with maximum tree heights of < 5 m (n = 381);
• NFI plots with stock per hectare = 0 (only calculated for trees with a diameter at breast height (dbh) of > 7 cm; n = 278).

In total, 1315 NFI plots were removed, equivalent to 3.68% of all plots. Subsequently, the NFI forest decision was compared to the forest cover map including 34,429 NFI plots (12,491 inside and 21,938 outside of the forest area).

2.4. Workflow

To calculate the forest cover, orthoimages and nDSMs as well as Sentinel-2 time series were used (see Figure 4). The challenge was to develop a method which could be flexibly adapted to the various aerial photographs and Sentinel-2 granules. Depending on the aerial image or satellite imagery, the spectral composition and lighting conditions varied.

To derive a vegetation mask from the orthoimages, 5000 pixels inside and 5000 pixels outside the forest area were selected from 200 randomly chosen orthoimages for each mission. Administrative data about forest ownership types served as a forest mask. Based on the specific reflective properties of vital vegetation, parameters can be derived that serve to differentiate between inhabited surfaces and living or dead vegetation. The best known of these parameters is the normalized difference vegetation index (NDVI), which is calculated from the ratio of the values of the red (R) and near infrared (NIR) spectral range [39]. Besides NDVI, greenness and brightness were calculated for each pixel (see Formulae (1)–(3)).

$$\mathrm{NDVI}= \frac{\mathrm{NIR} - \mathrm{Red}}{\mathrm{NIR} + \mathrm{Red}}$$

(1)

$$\mathrm{Greenness}= \frac{\mathrm{Green}}{\mathrm{Blue} + \mathrm{Red}}$$

(2)

$$\mathrm{Brightness}= \sqrt{\mathrm{Blue}² + \mathrm{Green}² + \mathrm{Red}²}$$

(3)

Since the absolute values for these parameters differed greatly between the flight missions, we used a Gaussian mixture model (GMM) to select the parameters corresponding to vegetation. The parametrization of a GMM is based on the expectation–maximization (EM) algorithm for fitting the Gaussian distributions [55]. The EM algorithm iteratively approximates the ideal model of Gaussian distributions underlying the observed data. For initializing the EM approach, we assumed five components as land cover classes to be present inside the forest areas of Baden-Württemberg: coniferous trees, broad-leaved trees, vegetated non-wood areas, non-wood areas without vegetation (i.e., bare soil, bedrock, paved forest roads) and shadow. The basic parameters of the GMM were the number of components (n = 5), covariance type (“full”) as the correlation of the used indices was in an acceptable range, maximum number of iterations (n = 100) and random state (n = 3) for reproducible results. The results of the GMM were 5 classes, which differed in terms of their value range of NDVI, greenness and brightness. The classes corresponding to vegetation were characterized by a high NDVI and greenness and were selected accordingly. The selection was visually verified by means of the orthoimages.

The second vegetation mask was derived from Sentinel-2 data. In order to eliminate all clouds and to obtain a homogeneous spectral composition for all granules, a time series was applied. From all Sentinel-2 scenes, the NDVI and subsequently the median NDVI of all scenes were calculated. The vegetation layer was created by setting a threshold value. This was possible due to the time series that achieved consistent NDVI values across Baden-Württemberg. By defining a threshold value (>0.7) for vegetation by visual interpretation of the orthoimages, a generalized vegetation layer of Baden-Württemberg was created covering forests along with shadow areas as well as meadows and green agricultural fields, without including settlement areas such as paved areas or rooftops. To generate a tree layer, the two vegetation layers were combined and compared with heights from the nDSM. Consequently, the tree layer consists of pixels, which were classified as vegetation based on orthoimages or Sentinel-2 data and have a minimum height of 3, 5 and 8 m determined by the nDSM. The forest cover map was calculated on the basis of this tree layer according to the forest definition by the German NFI, whereby the criteria canopy closure (>50%), minimum width (10 m), minimum size (0.1 ha) and minimum height (5 m) were considered.

The median of the NDVI time series of Sentinel-2 (Level2A) data was processed using the Google Earth Engine (GEE). The GEE computing platform offers a catalog of satellite imagery, geospatial datasets and planetary-scale analysis capabilities [56]. All further calculations were performed using a set of Python libraries for raster processing and analysis, GIS and scientific computing.

2.5. Comparison of Forest Cover Map and NFI Data

The forest cover map was compared to the forest decision (presence/absence of forest) of the NFI plot data (see Figure 5). To meet this objective, the confusion matrix, Kappa value and overall accuracy were calculated. NFI plot data served as reference data, whereas the forest cover map was used as a prediction. The threshold for forest presence was set to 50% forest cover, meaning that NFI plots with a forest cover of > 50% were classified as forest. Although the NFI plot data were filtered for this purpose, the remaining NFI plots do not completely correspond to forest cover. In the German NFI, minimum tree height is not defined, whereas it was an essential parameter for the derivation of the forest cover map. Furthermore, NFI data show land use, while the forest cover map represents land cover.

Based on the confusion matrix, differences between the forest cover map and NFI forest decision were identified and quantified. Therefore, the following data were used:

• Median nDSM heights within NFI plots: identification of forest plots < 5 m;
• Forest cover (%) of each plot: identification of plots at forest borders;
• Land use from the administrative layer: information about land use for each plot;
• Orthoimages: by visual interpretation, reasons for differences between the two datasets were identified and categorized.

3. Results

3.1. Forest Cover Map

As a result of this study, a wall-to-wall forest cover map for 99.71% of Baden-Württemberg was developed. It covers a total forest area of 13,222.1 km² with a spatial resolution of 1 m (see Figure 6). For 103.2 km² of the area (0.29% of Baden-Württemberg), no information (NoData) about forest cover could be generated due to missing aerial images, as described in Section 2.2.1. The final outcome of this study represents the first forest cover map of Baden-Württemberg in this very high level of detail. The distribution of forest cover appears to be reasonable. In particular, the high percentages of forest in the Odenwald, Black Forest and Swabian Alps as well as non-forested valleys in the Black Forest are clearly visible (see Figure 6). According to the forest cover map, approximately 13,222.1 km² (37.12%) of Baden-Württemberg is covered by forest. This is highly consistent with the results of the latest NFI [41], which reported a forest cover of approximately 13,319.6 km² (±157.8 km²), equal to 37.26% (±0.44%).

3.2. Comparison of NFI Plot Data and Forest Cover Map

The spatial agreement between the NFI data (based on a visual interpretation of forest presence/absence) and the forest cover map (based on remote sensing data) was 95.21% overall. The Kappa value showed an almost perfect agreement with a value of 0.90. The confusion matrix is shown in

Table 2. Confusion matrix of NFI plot data and forest cover map (fcm). OA = overall accuracy, UA = user’s accuracy, PA = producer’s accuracy, κ = Kappa value.

		NFI Plot Data
		Non-Forest	Forest	UA (%)
fcm	non-forest	21073	783	96.41
	forest	865	11708	93.12
	PA [%]	96.06	93.73
	OA [%]			95.21
	κ			0.90

. The producer’s and user’s accuracies for both forest and non-forest were very similar. Non-forest achieved, with 96.06–96.41%, a slightly higher producer’s and user’s accuracy than the forest class with a producer’s and user’s accuracy of 93.12–93.73%.

3.3. Analysis of the Differences between NFI Plot Data and Forest Cover Map

The forest cover map based on remote sensing and the NFI forest decision differed by 4.79%, which is 1648 NFI plots in total. Each NFI plot, where the NFI forest decision was not corresponding with the forest cover layer (= “incongruent NFI plot”), was analyzed on the basis of median nDSM height, forest cover percentage and land use data. In combination with visual interpretation, the following reasons for the incongruence between the two datasets could be categorized (see Figure 7):

• Edge effects: Incongruences occurring at forest borders. Plots at forest borders were determined to have a forest cover of 1–99%.
• Forest definition: The NFI forest definition does not include a minimum tree height. For this reason, during the NFI, all stocked areas were classified as forest, regardless of tree height. In contrast, the forest cover layer defines areas with a height of < 5 m as non-stocked.
• Land use: Actual land cover does not correspond to land use. There can be non-stocked forest areas (e.g., storm damages) or stocked non-forest areas (e.g., orchards).
• Layer errors: differences due to errors in the remote sensing-based forest cover layer (e.g., errors in image matching lead to wrong nDSM heights).

Based on the defined categories, each incongruent NFI plot was assigned to a specific category. The distribution of these NFI plots is illustrated in Figure 8. It can be observed that edge effects are evenly distributed all across the study site. The category forest definition occurs most often in the northern Black Forest which is characterized, partly, by a high proportion of openings in forests. In contrast to the category forest definition, the category land use is slightly stronger when represented outside of forest-dominated regions. Only three errors in the forest cover layer were identified.

Statistically, the categories edge effect and forest definition are the most common, with 2.08% and 2.04% of all NFI plots, respectively. The category land use represents 0.66% overall, whereas layer errors occur only in three cases, yielding a presence of 0.01%. To analyze the proportion of forest presence/absence within the categories, the NFI forest decision served as a reference. Forest presence/absence is distributed unevenly within the categories. Edge effects occurred in 1.85% of the NFI-based non-forest plots and only in 0.23% of the forest plots, whereas all NFI plots categorized to the class forest definition were situated in the forest. From the category land use, 0.65% of the NFI plots were located outside the forest. All errors in the forest cover map (n = 3) occurred outside the forest areas (see Figure 9). The NFI plots categorized as errors in the forest cover map are shown in Figure 10. Two of them are lying in agricultural fields, whereas one plot is located on a rooftop of a building.

4. Discussion

4.1. Forest Cover Map

In the present study, the first wall-to-wall forest map of Baden-Württemberg was produced. We agree with the statement of the authors of [1], who argued that the remotely sensed forest area mapping within the inventory cycles of an NFI provides information that extends the sample-based NFI information content and improves the reporting capacity. The derivation of the forest cover map is almost completely automated, meaning that the applied methods are suitable for area-wide forest mapping. Only for the creation of the vegetation layer is a selection of GMM classes representing vegetation necessary. Sentinel-2-derived time series together with aerial images enabled spatially and temporally explicit depictions of the forest area. By picking the median NDVI of the Sentinel-2 time series, clouds as well as cloud shadows could be removed. Variations in illumination as well as differences in the reflection of vegetation due to different acquisition times were minimized. The multisensory approach combined the advantages of both sensors: Sentinel-2 data performed better inside of forests due to less shadow effects compared to aerial images; orthoimages complemented the Sentinel-2 data by their high resolution. Thereby, single trees outside the forest area and very small forest patches could be detected. Similar to [7], the spatial resolution of the forest cover map is 1 m. As described in [30], spatially detailed data are preferred due to the following reasons: First, applying coarse-resolution data results in a local overestimation of the dominant land cover class and consequently in lower fragmentation rates [57]. For example, road corridors as an important contributor to forest fragmentation [58] can only be detected if the images have a sufficiently fine spatial resolution. Second, due to the average size of forest patches in Europe, an accurate analysis of forest patches requires data with a spatial resolution finer than 100 m [30]. Aerial images are affected by changing illumination conditions both between and within images [54]. This affects the reliability and robustness of spectral metrics derived from aerial images. We managed to compensate the effect of changing illumination as well as the influence of the different settings between the flight missions (see Appendix A) by introducing the GMM and selecting the resulting GMM classes corresponding to vegetation. With the developed methodology, it was possible to map the forest cover in a consistently high quality over the entire area of Baden-Württemberg. The results of this study show that a forest cover map can be derived by aerial images acquired in numerous photogrammetric image flights with varying times of flight and camera types as well as different flight and camera settings. Across the forest cover map, errors are equally distributed, indicating that they are unaffected by the differences between the aerial image flights or the Sentinel-2 granules.

As the methodology is entirely based on the analysis of the spectral properties of land cover, the resulting map represents a forest cover map rather than a forest (land use) map. Furthermore, the combination of spectral and height information basically only captures the visible canopy. Objects below cannot be considered (e.g., parking lots or camping sites in the forest). Furthermore, the forest cover map contains also, for example, field shrubs, hedges, parks, orchards, short-rotation plantations and energy wood areas. On the other hand, temporary or permanently unstocked areas (e.g., timber stockyards, forest roads), which by definition are forest, are not recorded. Depending on the applicable forest definition and task, these areas must be added or subtracted [39].

4.2. Comparison NFI Plot Data and Forest Cover Map

With a few exceptions, the present forest cover map agrees well with the NFI data. The high overall accuracy of 95.21% of the forest cover layer indicates that the applied methods and results are promising. However, it should be considered that this kind of accuracy assessment takes into account all plots including plots at forest edges, plots with a canopy height of <5 m and plots where land cover does not match the actual land use. Since the forest cover layer only has the claim to cover forest without considering land use, the accuracy of the forest cover layer is much higher than the calculated overall accuracy. This is also evident in the small number of incongruences caused by layer errors. As mentioned in [7], a comparison with other studies is difficult. Studies, e.g., mapped forest with different input data types, had substantially smaller study areas or used varying forest definitions. To compare the results with other studies on forest cover, the studies mentioned in the Introduction can be reconsidered. In the following, a brief comparison is drawn with studies that also used NFI data as a reference: Ref. [7] obtained, with DSMs from image-based point clouds, an overall accuracy of 97% for the forest cover map of Switzerland taking the land use criterion into account; Ref. [18] used Sentinel-2 data for tree cover mapping for two study sites in Europe and achieved an overall accuracy of up to 90%; Ref. [1] used Landsat data to map the forests of Canada and achieved an overall spatial agreement of 84%.

Mapping of forest areas is affected by the categorical, spatial and temporal elements involved [1]. Different forest definitions as well as data processing and analysis methods lead to different mapping results [59]. In our study, we identified four reasons for incongruences between the NFI forest decision and forest cover layer: (a) edge effects, (b) different forest definitions, (c) land cover does not match land use and (d) layer errors.

The authors of [7] identified forest borders as one of the main reasons for mismatches between the remotly sensed forest cover map and NFI data. They revealed that in densely forested or entirely non-forested areas, the forest cover map has a higher accuracy than at forest borders. Similar to the aforementioned study, forest borders explained 2.08% of all incongruences. If an NFI tract corner was placed onto a road, no inventory was performed. Therefore, the category edge effect includes not only NFI plots on forest borders but also NFI plots on roads intersecting forests. Furthermore, the co-registration of field data with remote sensing data causes information uncertainties due to position inaccuracies. The NFI plots have an unknown positioning error due to inaccuracies of the GNSS signals below the forest canopy. This can lead to a shift between the observations in the field and in the remote sensing data [46]. Consequently, the category edge effects may also contain NFI plots with position errors.

In addition, the NFI does not define minimum tree height, which led to incongruences between the NFI forest decision and forest cover map. Forest definition is basically a subcategory of the category land use. The categories might overlap, since the assignment of non-stocked areas such as forest gaps to the category forest definition is based on the fact that land use is forest. If tree height is too low, but the plot is categorized by the NFI as forest, land use is forest, but land cover is not. That is the reason why these incongruent plots are all inside the forest. Based on the nDSM, a minimum vegetation height can be identified, which was applied in the form of a threshold value eliminating all objects that lie below this minimum height. When setting such a threshold value, it should be noted that the lower it is set, the more difficult it is to separate forest vegetation from other forms of vegetation. If the minimum height is set low, the likelihood of confusing forest vegetation with high-growing herbaceous vegetation, e.g., corn fields, large grasses and tall perennials, increases. Small trees in plantations and natural regeneration are not shown in the forest cover map due to their small size [39].

Generally, leisure facilities such as parks and sports fields, infrastructure such as roads (avenues) and cemeteries and agricultural areas such as green strips or orchard meadows are not considered as forests. Furthermore, in the category land use, there are 53 NFI plots where the forest decision is not verifiable. This might be explained by errors during preclarification. New, previous non-forest tracts, which are located with at least one corner in the forest, were overseen.

In addition, the defined minimum forest area of 0.1 ha takes also very small forests into account. Most of these areas are not determined as forest in terms of land use, but in the forest cover map, they are declared to be stocked.

The results show that there are almost no errors in the forest cover layer (n = 3). If the land surface depicted in aerial images is too homogeneous, image matching may be incorrect, resulting in artifacts in the orthoimages and wrong heights in the elevation models. Depending on the selected class of the Gaussian mixture model, it is possible that green, gray or light roofs in the orthoimages will be classified as vegetation. Most of the invalid areas are discarded as soon as the minimum height is considered using the nDSM. However, in the case of high-non-vegetation objects such as houses or high-voltage lines, non-vegetation may be included in the forest layer. Due to the Sentinel-2 pixel size of 10 × 10 m, parts of non-vegetation areas such as corners of rooftops can be included in the forest cover map.

The comparison of the remote sensing-based forest cover and the percentage forest cover of Baden-Württemberg reported by the latest NFI [41] shows a high correspondence with a slight difference of 0.14 percent points. Mismatches can be explained with the categories analyzed in Section 3.3. The fact that the reported forest area is higher than the calculated forest cover leads to the assumption that, in Baden-Württemberg, there are more forest areas with tree heights of <5 m than stocked areas outside of forests. This might be explained partially by the high proportion of openings in the northern Black Forest (see Figure 8).

5. Conclusions

Our study developed the first forest cover map of Baden-Württemberg using remote sensing data. The comparison with the NFI data showed that the applied methods and results are promising. Satellite-derived time series in combination with aerial images enable spatially and temporally explicit depictions of stocked areas in a high level of detail. Due to the free availability of the Sentinel-2 products as well as the regular image flights of the LGL, data availability is guaranteed on a long-term basis. The tree and forest cover maps are currently available for the years from 2011 to 2019 and will be updated regularly. In the near future, this will enable the elaboration of time series analysis and thus will allow long-term and area-wide monitoring of tree and forest cover over Baden-Württemberg. Not only in the future but already now, forest cover maps provide valuable information for multiple research questions in biodiversity studies, forest resource analysis, climate change research and other forest-related studies.