The Efficacy Analysis of Determining the Wooded and Shrubbed Area Based on Archival Aerial Imagery Using Texture Analysis

Abstract

Open areas, along with their non-forest vegetation, are often threatened by secondary succession, which causes deterioration of biodiversity and the habitat’s conservation status. The knowledge about characteristics and dynamics of the secondary succession process is very important in the context of management and proper planning of active protection of the Natura 2000 habitats. This paper presents research on the evaluation of the possibility of using selected methods of textural analysis to determine the spatial extent of trees and shrubs based on archival aerial photographs, and consequently on the investigation of the secondary succession process. The research was carried out on imagery from six different dates, from 1971 to 2015. The images differed from each other in spectral resolution (panchromatic, in natural colors, color infrared), in original spatial resolution, as well as in radiometric quality. Two methods of textural analysis were chosen for the analysis: Gray level co-occurrence matrix (GLCM) and granulometric analysis, in a number of variants, depending on the selected parameters of these transformations. The choice of methods has been challenged by their reliability and ease of implementation in practice. The accuracy assessment was carried out using the results of visual photo interpretation of orthophotomaps from particular years as reference data. As a result of the conducted analyses, significant efficacy of the analyzed methods has been proved, with granulometric analysis as the method of generally better suitability and greater stability. The obtained results show the impact of individual image features on the classification efficiency. They also show greater stability and reliability of texture analysis based on granulometric/morphological operations.

1. Introduction

In Poland, since the 1990s, a gradual cessation of agricultural use, hence the secondary succession process, has been observed in some areas. This is also observed in valuable natural areas, including Natura 2000 habitats, and results in deterioration of biodiversity and the habitat’s conservation status. This process can have different characteristics and dynamics in different areas [1]. The course of succession is influenced by, e.g., the state of abandoned land (e.g., black fallow, mowed meadow, stubble), soil fertility in nutrients, recently grown plants, the location of the ground (including proximity to the forest), slope, and insolation [2,3,4]. This knowledge plays an important role in proper planning of active protection of the Natura 2000 habitats.

The secondary succession process can now be effectively monitored using LiDAR (Light Detection and Ranging) data [5,6] and hyperspectral imagery (detection of succession species) [7]. However, learning about the dynamics of this phenomenon in the past (e.g., in the 1970s, 1980s, and 1990s of the 20th century) is only possible using archival materials—mainly aerial photographs, which are often the only reliable source of information about the land cover. The most common method of obtaining ranges of trees and shrubs based on archival aerial photographs is photo interpretation [8,9,10,11,12,13]. Sometimes, stereodigitalization is also used [14,15,16,17,18]. However, these methods are labor-intensive and time-consuming, which is why methods of automation of this process are sought. Sometimes, the classic spectral classification of the image is used [19], but in the case of archival images (often black and white), it is not always possible. Another possibility for the study of the succession process of trees and shrubs is to use the technique of dense image matching (DIM) [20]. However, according to current research in this area [20], for the DIM technique to yield good quality results, archival data must meet several conditions: (1) The scale of images should not be less than 1:13 000 for analog photos, or GSD (Ground Sample Distance) should be smaller than 25 cm for digital photos; (2) the date of obtaining the photos should provide a picture of vegetation in full development; and (3) photos should have good radiometric quality, i.e., small image blur, adequate contrast, and no mechanical damage (in the case of analogue photos).

A new approach to the study of the dynamics of the succession process based on archival aerial photographs may be the application of textural analysis of the image.

Texture is one of the most important spatial features of selected terrain coverage classes. Since it does not have an unambiguous mathematical definition, in practice, image processing uses a variety of ad hoc texture analysis methods. The gray level co-occurrence matrix (GLCM) [21,22,23], fractal analysis [24], discrete wavelet transformation [25], Laplace filters [26,27,28], Markov Random Fields [29,30], and granulometric analysis [31,32] are several such methods that should be mentioned here. Also, promising attempts have been made to use convolutional neural networks to take spatial features of images into account [33,34].

Wooded and bushy areas are a class with a distinct texture. This is mainly due to the shape of the tree crowns and also to their height. In areas with dense trees/bush coverage, typical grainy texture results from the alternate occurrence of better and worse illuminated crown fragments (due to the orientation of the surface of the tree crown relative to the light source). In areas with less dense coverage, the shadow image of individual trees or shrubs also plays an important role.

Texture is especially important for images of very high resolution. As the size of the pixel (GSD) increases, the clarity of texture decreases, and at the same time, its significance as an interpretive/classification feature does too [35,36].

The texture, although potentially significant in the classification of satellite image content, used alone might not give satisfactory results, because different LULC (Land Use/Land Cover) classes may have similar textures [37]. Therefore, the use of textural features as a complement to spectral features gives a much better result—in the spectro-textural approach. This is indicated by examples of models using different texture analysis methods in combination with spectral data [22,23,24,33,34,35,36,38,39,40,41,42,43].

Most examples of the spectro-textural approach are based on multispectral data, including near-infrared channels (often also mid-infrared), which are crucial for distinguishing vegetation from other terrain coverage classes, as well as for differentiating types of vegetation. In the case of archival aerial photographs, these type of data are often not available. Evidently, the repository of multispectral satellite data dates back to the 1970s, but these are medium-resolution images, unsuitable for the analysis of textures of wooded areas or for large-scale studies. On the other hand, aerial photographs usually have adequate spatial resolution, while the registered spectral band images are less favorable for texture analysis. Most often, these are images in natural colors or panchromatic (black and white) photographs.

The main objective of the presented research was to analyze the possibility of using archival aerial imagery to determine the wooded and shrubbed areas, using selected methods of texture analysis, and with limited spectral data. Two methods of texture analysis were selected for the study, the effectiveness of which has been proved by previous studies: GLCM [21,23] and granulometric analysis [31,32,33,36]. In the literature can be found many studies on the use of these methods in the analysis of satellite and aerial images, but there is almost no research on archival aerial photos.

2. Brief Presentation of Selected Methods of Texture Analysis

Here follows a brief presentation of two texture analysis methods used for this research: GLCM and granulometric analysis. Both methods can be used to analyze the whole image, as well as for local analysis for individual pixels. In the first case, we obtain a single value (or a set of values) describing the feature (or features) of the whole image. This may, for example, determine what the image presents (e.g., the dominant type of coverage/land use) [44,45]. The local analysis consists of analyzing a subset of the image defined by a determined neighborhood of a single pixel (e.g., all pixels in the assumed radius). Values obtained in this way are attributed to individual central pixels. This indicates the feature of texture in the environment surrounding each central pixel (defined by the range of the analyzed neighborhood). In practice, this consists of creating new images, where pixels are given values depending on the individual features obtained in the analysis.

2.1. GLCM—Grey Level Co-occurrence Matrix

The GLCM method of texture analysis was first introduced by Julesz [21]. It consists of counting the pairs of neighboring values present in the image (or in its subset) and comparing them in the (grey level co-occurrence) matrix. For the matrix thus created, selected indicators are calculated to analyze various aspects of the texture of the image (or its subset). Works of Haralick et al. [23], in which numerous indicators were presented, played an important role in developing and popularizing this method. Therefore, GLCM indicators are frequently referred to as Haralick indicators (or features). This method is highly effective [46], which has been demonstrated in many publications [47,48], although different authors suggest using different sets of Haralick features [38,39,41]. One of the main advantages of GLCM is the large number of indicators describing various aspects of the texture, enabling to analyze it from different angles. In these studies, we used a set of eight basic Haralick indicators (formulas according to the OTB CookBook [49]).

$$\mathrm{Energy}={\displaystyle \sum }_{\mathrm{i},\mathrm{j}}\mathrm{g}{\left(\mathrm{i},\mathrm{j}\right)}^2$$

(1)

Energy is a description of the uniformity of texture.

$$Entropy={{\displaystyle \sum }}_{i,j}g\left(i,j\right)lo{g}_{2}g\left(i,j\right),or0ifg\left(i,j\right)=0,$$

(2)

Entropy is a measure of randomness of pixel values distribution. It is the highest when all the frequencies g(i,j) are equal, and smaller when they are unequal. Heterogeneous images (subsets of image) will result in a higher entropy value.

$$Correlation={{\displaystyle \sum }}_{i,j}\frac{\left(i-\mu \right)\left(j-\mu \right)g\left(i,j\right)}{{\sigma }^2}$$

(3)

Correlation describes how a pixel is correlated to its neighborhood. Regions of similar gray level will result in higher correlation values.

$$InverseDifferenceMoment={{\displaystyle \sum }}_{i,j}\frac{1}{1+{\left(i-j\right)}^2}g\left(i,j\right),$$

(4)

Inverse difference moment (IDM) measures a homogeneity of an image (its subset). Homogeneous images will result in high IDM values.

$$Inertia={{\displaystyle \sum }}_{i,j}{\left(i,j\right)}^{2}g\left(i,j\right),$$

(5)

Inertia (or contrast) measures a contrast of pixel values between a given pixel and its neighborhood. Images with a large amount of variation will result in high inertia.

$$ClusterShade={{\displaystyle \sum }}_{i,j}{\left(\left(i-\mu \right)+\left(j-\mu \right)\right)}^{3}g\left(i,j\right),$$

(6)

Cluster shade is a measure of the skewness of the matrix. It describes the perceptual concepts of uniformity [50].

$$ClusterProminence={{\displaystyle \sum }}_{i,j}{\left(\left(i-\mu \right)+\left(j-\mu \right)\right)}^{4}g\left(i,j\right),$$

(7)

Cluster prominence is a measure of asymmetry. The less symmetric regions will result in higher cluster prominence values.

$$Haralic{k}^{\prime }sCorrelation=\frac{{{\displaystyle \sum }}_{i,j}\left(i,j\right)g\left(i,j\right)-{\mu }_t^2}{{\sigma }_t^2},$$

(8)

where $\left(i,j\right)$ is the matrix cell index, $g\left(i,j\right)$ is the frequency value of the pair having index $\left(i,j\right)$, $\mu ={{\displaystyle \sum }}_{i,j}i\ast g\left(i,j\right)={{\displaystyle \sum }}_{i,j}j\ast g\left(i,j\right)$ (due to matrix symmetry) and means weighted pixel average, $\sigma ={{\displaystyle \sum }}_{i,j}{\left(i-\mu \right)}^2\ast g\left(i,j\right)={{\displaystyle \sum }}_{i,j}{\left(j-\mu \right)}^2\ast g\left(i,j\right)$ (due to matrix symmetry) and means weighted pixel variance, and μ_t and σ_t are the mean and standard deviation of the row (or column, due to symmetry) sums.

2.2. Granulometric Analysis

The granulometric analysis in a classical form consists of sequentially executing a series of morphological openings of a binary image using structuring elements of sequentially increasing size, and then in the calculation of derivative images demonstrating the differences [31]. By specifying the size of each differential image (defined as the sum of all pixel values in the image) [51], we are able to determine the number of texture grains of specific sizes found in the source image [31]. Due to the nature of the opening operation (removal of relatively small image elements with values greater than the background value), this applies to bright texture elements. Analogically to the analysis based on the sequence of morphological openings, an analysis based on morphological closing operations can be performed. Due to the dual nature of the closing operation (in opposition to the opening operation), we obtain complementary information about the dark texture elements. This type of analysis is sometimes referred to as anti-granulometric [35]. However, in this article, both versions of the analysis (opening- and closing-based) are referred to as granulometric analysis. This type of global analysis assigns a set of features called the granulometric density function to the whole image.

The use of local granulometry is presented in [32,52]. It consists of assigning the granulometric density function to individual pixels on the basis of the analysis of their neighborhood. In practice, the effect of local granulometric analysis is usually a set of images in which pixel values represent successive values of the granulometric density function. Such images are called granulometric maps [30,50]. The usefulness of granulometric maps in the LULC classification has been proven, e.g., in [6,35,40,51]. It is worth noting that this type of granulometric analysis in many respects resembles a morphological profile [53,54]; however, it differs significantly in other respects [35,43].

Among the important advantages of granulometric analysis, a multiscale should be indicated: As a result of successive morphological operations with an increasing size of the structuring element, subsequent values (granulometric maps) present information about the presence of texture grains of different sizes.

Another important feature of this method is resistance to the so-called edge effect [55]. Most texture analysis methods are based on the analysis of the spatial frequency of the image, as a result of which the edges of the objects, even those with a low texture, are marked with values indicating high texture. This may lead to a reduced accuracy of the classification [35,36,55]. The granulometric analysis is based on a different principle: It removes whole elements of the image and analyzes their number; as a result of this, it does not show the edges of objects as areas with high texture.

3. Study Area

The study area is located in the south of Poland, in the Silesian Voivodeship, near Czestochowa city (50°45″ N; 19°17″ E) (Figure 1). This area is a part of the Olsztyn-Mirów wildlife refuge, a Natura 2000 protected area (PLH240015), established by the decision of the European Commission of 12 December 2008 No. 2009/93/EC (OJ EU L 43 of 13.02.2009) [56].

This area is characterized by large habitat diversity. Non-forest habitats present there are of particular importance. They are associated with limestone rocks, with numerous rare and endangered thermophilic species of plants and invertebrates (including the scarce large blue (Phengaris teleius), a species from Annex II of the Habitats Directive). A number of species reaches here the northern edge of their geographic range. In total, there are 13 types of habitats identified in the refuge with Annex I of the Habitats Directive (including 4 priority habitats) and 12 species of plants and animals from Annex II of the Habitats Directive [57].

The threat to the non-forest habitats of this area is primarily the process of secondary succession, which is the result of the abandonment of pastoralism, the lack of grazing and mowing of meadows, which has been progressing gradually since the 1990s. In order to protect valuable habitats at the Olsztyn-Mirów wildlife refuge within the framework of the LIFE11 NAT/PL/432 project “Protection of valuable non-forest habitats of the Orle Gniazda Landscape Park” (2012–2016), a number of active protection measures, such as the planting of bushes and trees and controlled sheep grazing, have been undertaken [58].

The part of the area selected for analysis (25 km²) covers a territory which is varied in terms of terrain relief and land use, and of the nature of the secondary succession process. It allows to identify possible problems in the automation of the process of determining the succeeding ranges of trees and shrubs using textural analysis methods.

4. Materials and Methods

The main stages of the development methodology included:

• Analysis and selection of archival aerial imagery;
• Development of orthophotomaps;
• Creation of simulated panchromatic (P) images;
• Determination of the wooded and shrubbed areas using textural analysis, including:
  • Selection of methods and parameters of texture analysis;
  • Creation of spectro-textural data sets;
  • Classification using support vector machine (SVM); and
  • Additional processing.
• Preparation of reference data and evaluation of the accuracy of determining the wooded and shrubbed areas; and
• Analysis and comparison of results.

The methodology scheme is presented in Figure 2.

4.1. Data

The aerial photographs for the study area were obtained from the State Geodetic and Cartographic Resource in Poland. They were photos from 1971, 1996, 2003, 2009, 2012, and 2015. They were both analogue and digital, black and white (panchromatic—P) and color (RGB). They differed in parameters, scale (from 1:13 000 to 1:26 000) or GSD (24 and 25 cm), and the acquisition date (different phenological periods). The basic parameters characterizing these photographs are presented in

Table 1. The characteristics of aerial photographs used for analysis.

Date	Number of Photos	GSD or Scale	Camera	Focal Length	GPS/INS	Aerotrian-Gulation (EO)	Type
11.08.1971	12	1:18 000	RC51	210.20 mm	NO	NO	P
30.05.1996	4	1:26 000	RC20	152.97 mm	NO	NO	RGB
24.05.2003	14	1:13 000	LMK	152.30 mm	NO	NO	P
26.04.2009/ 29.04.2009	14	1:14 000	RC30	153.81 mm	YES	YES	RGB
25.03.2012	10	24 cm	UltraCamXp	100.50 mm	NO	NO	RGB
08.08.2015	10	25 cm	UltraCamXp	100.50 mm	NO	NO	RGB, CIR

, and in Figure 3, selected parts of photos are visualized.

Based on the obtained images, orthophotomaps with GSD = 0.5 m were prepared. In the case of scanned images, first they were pre-processed in SAPC software according to the methodology developed by Salach [59]. This methodology transforms scans to a form similar to the images captured by a digital camera. Next, approximate values of exterior orientation (EO) parameters were obtained using Agisoft Photoscan software. Then, these values were used as approximate EO values in Trimble INPHO, which allowed the process of image block bundle adjustment to occur. A more detailed description of the process of developing archival materials is in [20].

4.2. Simulation of Panchromatic Images

Due to the importance of spectral data for spectro-textural classification, one of the aspects of this study was to compare the effectiveness of classification depending on the type of images available. However, the test pictures, apart from the type, also differ in the radiometric quality, which can undoubtedly be relevant to the accuracy of the quality of texture analysis and, as a result, to the classification performed. In order to estimate the importance of the type of photo regardless of the quality of the material itself, panchromatic images were simulated for RGB images.

Aerial panchromatic photographs images were most often made with the use of a yellow (minus-blue) filter blocking blue radiation [60,61]; so only green and red radiation was registered in panchromatic images. Therefore, simulated panchromatic images ($P_s$) were calculated by averaging the values of red ($R$) and green ($G$) bands:

$$P_s=\frac{R+G}{2},$$

(9)

Simulated panchromatic images were generated for all RGB images—that is, for data from 1996, 2009, 2012, and 2015.

4.3. The Determination of Wooded and Shrubbed Areas

The determination of wooded and shrubbed areas was based on spectro-textural classification. It was made using the support vector machine (SVM) method on combined spectral data (original CIR, RGB, or P and simulated P) and texture data (Haralick indicators or granulometric maps).

4.3.1. Selection of Source Images for Texture Processing

Single grayscale images are proposed for texture processing. In the case of P images (both original and simulated), the choice was obvious—to choose these P images. In the case of multispectral images (CIR and RGB), the image of the first principal component was selected for texture processing. This image by definition represents the greatest variation of the data set. Previous research also shows that texture analysis based on the image of the first principal component gives the best distinction of the texture of selected LULC classes—in addition to the near-infrared image [62], unavailable for most of test images in this study.

4.3.2. Texture Analysis

Texture analysis was done independently using two methods: GLCM and granulometric analysis. The two proposed methods were chosen for the following reasons. GLCM is one of the most popular methods of texture analysis, according to the literature analysis still very often used in the analysis of remote sensing data [22,23,41,47]. However, granulometric analysis is a relatively well-known method of texture analysis, but its use in remote sensing is negligible. Earlier studies show that both methods are highly effective, also when compared to other texture descriptors. In addition, these methods are highly reliable and are easy to implement. In both cases, it was a local analysis performed in three variants, depending on the size of the analyzed neighborhood:

• 10-pixel radius;
• 15-pixel radius; and
• 20-pixel radius.

GLCM analysis was performed using eight Haralick features, according to the formulas given in Equations (1)–(8). As a result, eight images were obtained, which were treated as complementary texture information, basing on all four displacement vectors. The OTB Monteverdi software was used for the calculations.

The granulometric analysis based on both opening and closing operations. For each image, 6 granulometric steps (opening or closing) were made, resulting in 12 granulometric maps (6 for opening and 6 for closing), representing information about texture corresponding to the size of structuring elements (circle-shaped, from a 1-pixel to 6-pixel radius). In order to investigate which granulometric maps carry information of significance from the viewpoint of the analysis of the wooded and shrubbed areas, data sets with fewer numbers of granulometric maps were also prepared. As a result, data sets with the following numbers of granulometric maps were analyzed:

• From 1 to 3 for opening and closing (6 granulometric maps);
• From 1 to 4 for opening and closing (8 granulometric maps);
• From 1 to 5 for opening and closing (10 particle size maps); and
• From 1 to 6 for opening and closing (12 granulometric maps).

The calculations were carried out using BlueNote software [63].

4.3.3. Tested Variants

Variants (165) of data sets for classification were prepared, differing in at least one of the following features: Date of acquisition, type of image, texture analysis method, local texture analysis neighborhood size, or the number of granulometric maps used (in the case of granulometric analysis). These variants are summarized in

Table 2. Tested variants of data sets.

Year	Image Type	GLCM			Granulometric Analysis
		10	15	20	10	15	20
1971	P	p-1971-glcm-10	p-1971-glcm-15	p-1971-glcm-20	p-1971-x*-10	p-1971-x*-15	p-1971- *-20
1996	P	p-1996-glcm-10	p-1996-glcm-15	p-1996-glcm-20	p-1996-x*-10	p-1996-x*-15	p-1996- *-20
	RGB	c-1996-glcm-10	c-1996-glcm-15	c-1996 glcm-20	c-1996-x*-10	c-1996-x*-15	c-1996- *-20
2003	P	p-2003-glcm-10	p-2003-glcm-15	p-2003-glcm-20	p-2003-x*-10	p-2003-x*-15	p-2003- *-20
2009	P	p-2009-glcm-10	p-2009 glcm-15	p-2009-glcm-20	p-2009-x*-10	p-2009-x*-15	p-2009- *-20
	RGB	c-2009-glcm-10	c-2009-glcm-15	c-2009-glcm-20	c-2009-x*-10	c-2009-x*-15	c-2009- *-20
2012	P	p-2012-glcm-10	p-2012-glcm-15	p-2012-glcm-20	p-2012-x*-10	p-2012-x*-15	p-2012- *-20
	RGB	c-2012-glcm-10	c-2012-glcm-15	c-2012-glcm-20	c-2012-x*-10	c-2012-x*-15	c-2012- *-20
2015	P	p-2015-glcm-10	p-2015-glcm-15	p-2015-glcm-20	p-2015-x*-10	p-2015-x*-15	p-2015- *-20
	RGB	c-2015-glcm-10	c-2015-glcm-15	c-2015-glcm-20	c-2015-x*-10	c-2015-x*-15	c-2015- *-20
	CIR	cir-2015-glcm-10	cir-2015-glcm-15	cir-2015-glcm-20	cir-2015-x*-10	cir-2015-x*-15	cir-2015-x*-20

* x represents the number of subsequent granulometric maps (after opening and closing) used to create a given variant of the data set (from 3 to 6).

.

4.3.4. Performing the Classification

In total, a classification was made on 165 variants of data sets. Fifteen variants were tested for each of the source images: 3 variants based on GLCM analysis and 12 variants based on granulometric analysis (the larger number is due to testing 4 sub-variants depending on the number of granulometric maps used).

The classification was made using the support vector machine (SVM) using ArcGIS software. This method was chosen because of its high efficiency compared to other machine learning methods, such as random forest or maximum likelihood [64].

The training fields were developed as a result of manual digitalization based on the visual interpretation of the original spectral images. All variants developed within one date of acquisition were classified based on the same set of training data. The final effect of each classification was the binary mask of wooded and shrubbery areas, created as a result of aggregation of individual classes based on individual training fields. The basic parameters related to the implementation of the classification are presented in

Table 3. Basic parameters of the classification of test images.

Date	Number of Variants Tested	Number of Training Fields: Wooded and Shrubbed Areas/Other Classes (Number of Pixels in Brackets)
1971	15	6/12 (1812/12180)
1996	30	4/15 (3832/4825)
2003	15	6/15 (2653/9057)
2009	30	13/32 (1379/6032)
2012	30	7/18 (8531/9520)
2015	45	11/17 (4590/5445)

.

4.3.5. Additional Processing

After obtaining the binary masks of wooded areas, they were additionally processed using morphological closing by reconstruction [65]. This operation was aimed at increasing the accuracy by removing any discontinuities in the mask of the wooded areas. Such discontinuities resulted mainly from the occurrence of shaded areas, i.e., those whose spectral features could significantly differ from those of other parts of wooded areas. Direct allocation of these areas to the wooded area mask could lead to an over-classification of this class in other parts of the image, and thus to a significant increase in the excess error and, as a result, to a decrease in the accuracy of the classification.

Closing by reconstruction removes (closes) relatively small discontinuities (of a size smaller than the size of the structuring element) in binary masks; at the same time, it does not change the shape of other objects, unlike other filters of a similar type (simple closure, median filter, majority, etc.).

The operation used closing by reconstruction using a 10-pixel radius structuring element. It was applied to all 165 classification variants obtained and the accuracy of newly obtained masks was independently assessed. The processing was carried out in the BlueNote software [63].

4.4. Accuracy Assessment

Evaluation of the accuracy of determining the wooded and shrubbed areas using various approaches to textural analysis was carried out, comparing their results with the reference data, which were the results of detailed visual interpretation of orthophotomaps prepared on the basis of the same archival materials. The entire analysis area was digitized for accuracy assessment purposes. Orthophotomaps were developed with GSD = 0.5 m. The smallest separation was of 2 m²; that is, it covered about 3 × 3 pixels. Thanks to such detailed visual interpretation, it was possible—in the later stages—to indicate: (1) Which of the textural analysis variants gives the best results, and (2) what size of wooded and shrubbed areas can be automatically determined on the basis of individual archival materials.

The following accuracy indicators were used to assess the efficacy of individual approaches: Overall accuracy (OA), recall, precision, Cohen’s kappa coefficient, and F1-score [66,67]. In addition, the sum of the area resulting from errors of omission and commission (EO + EC) was analyzed to indicate the variant with the smallest error area.

5. Results

The results of the analysis of the accuracy of determining the range of trees and shrubs using spectro-textural classification are presented in

Table 4. Results of the accuracy assessment obtained for individual images and analysis variants (the best results are marked in green).

Variant	P						RBG (C)						CIR
	K	OA	F1	PA	UA		K	OA	F1	PA	UA		K	OA	F1	PA	UA
	1971
3-10	0.722	0.932	0.762	0.714	0.816		-	-	-	-	-		-	-	-	-	-
3-15	0.633	0.910	0.685	0.645	0.731		-	-	-	-	-		-	-	-	-	-
3-20	0.634	0.913	0.684	0.663	0.707		-	-	-	-	-		-	-	-	-	-
4-10	0.726	0.936	0.763	0.752	0.775		-	-	-	-	-		-	-	-	-	-
4-15	0.644	0.916	0.693	0.677	0.709		-	-	-	-	-		-	-	-	-	-
4-20	0.624	0.912	0.675	0.667	0.683		-	-	-	-	-		-	-	-	-	-
5-10	0.714	0.933	0.752	0.746	0.759		-	-	-	-	-		-	-	-	-	-
5-15	0.669	0.921	0.715	0.694	0.737		-	-	-	-	-		-	-	-	-	-
5-20	0.614	0.914	0.663	0.696	0.633		-	-	-	-	-		-	-	-	-	-
6-10	0.692	0.929	0.733	0.738	0.727		-	-	-	-	-		-	-	-	-	-
6-15	0.641	0.921	0.687	0.726	0.652		-	-	-	-	-		-	-	-	-	-
6-20	0.563	0.904	0.618	0.660	0.581		-	-	-	-	-		-	-	-	-	-
glcm-10	0.775	0.945	0.807	0.754	0.868		-	-	-	-	-		-	-	-	-	-
glcm-15	0.676	0.928	0.717	0.758	0.681		-	-	-	-	-		-	-	-	-	-
glcm-20	0.682	0.932	0.720	0.791	0.661		-	-	-	-	-		-	-	-	-	-
	1996
3-10	0.724	0.895	0.794	0.720	0.883		0.753	0.904	0.817	0.727	0.931		-	-	-	-	-
3-15	0.727	0.891	0.800	0.691	0.948		0.749	0.901	0.815	0.710	0.956		-	-	-	-	-
3-20	0.621	0.850	0.720	0.630	0.840		0.690	0.877	0.772	0.671	0.907		-	-	-	-	-
4-10	0.714	0.892	0.785	0.721	0.863		0.744	0.905	0.806	0.758	0.860		-	-	-	-	-
4-15	0.671	0.863	0.761	0.633	0.955		0.735	0.894	0.806	0.693	0.963		-	-	-	-	-
4-20	0.647	0.857	0.742	0.632	0.901		0.702	0.883	0.780	0.684	0.906		-	-	-	-	-
5-10	0.703	0.884	0.780	0.691	0.894		0.733	0.893	0.804	0.695	0.952		-	-	-	-	-
5-15	0.672	0.863	0.762	0.633	0.958		0.744	0.898	0.812	0.702	0.963		-	-	-	-	-
5-20	0.668	0.866	0.757	0.647	0.912		0.668	0.869	0.755	0.663	0.876		-	-	-	-	-
6-10	0.695	0.878	0.776	0.669	0.923		0.723	0.889	0.797	0.685	0.951		-	-	-	-	-
6-15	0.650	0.854	0.747	0.619	0.941		0.707	0.880	0.786	0.665	0.962		-	-	-	-	-
6-20	0.611	0.846	0.713	0.622	0.835		0.642	0.860	0.735	0.648	0.850		-	-	-	-	-
glcm-10	0.527	0.796	0.661	0.533	0.869		0.597	0.820	0.714	0.561	0.982		-	-	-	-	-
glcm-15	0.337	0.723	0.521	0.431	0.657		0.572	0.811	0.696	0.551	0.944		-	-	-	-	-
glcm-20	0.243	0.687	0.450	0.377	0.558		0.305	0.716	0.493	0.418	0.602		-	-	-	-	-
	2003
3-10	0.761	0.889	0.846	0.757	0.957		-	-	-	-	-		-	-	-	-	-
3-15	0.750	0.886	0.836	0.768	0.919		-	-	-	-	-		-	-	-	-	-
3-20	0.756	0.890	0.838	0.790	0.892		-	-	-	-	-		-	-	-	-	-
4-10	0.755	0.887	0.841	0.763	0.937		-	-	-	-	-		-	-	-	-	-
4-15	0.760	0.892	0.842	0.787	0.904		-	-	-	-	-		-	-	-	-	-
4-20	0.734	0.883	0.820	0.804	0.837		-	-	-	-	-		-	-	-	-	-
5-10	0.754	0.886	0.841	0.754	0.951		-	-	-	-	-		-	-	-	-	-
5-15	0.761	0.892	0.843	0.782	0.915		-	-	-	-	-		-	-	-	-	-
5-20	0.739	0.885	0.824	0.802	0.848		-	-	-	-	-		-	-	-	-	-
6-10	0.735	0.879	0.826	0.762	0.902		-	-	-	-	-		-	-	-	-	-
6-15	0.761	0.892	0.843	0.784	0.910		-	-	-	-	-		-	-	-	-	-
6-20	0.727	0.881	0.815	0.809	0.821		-	-	-	-	-		-	-	-	-	-
glcm-10	0.749	0.884	0.838	0.751	0.947		-	-	-	-	-		-	-	-	-	-
glcm-15	0.748	0.885	0.835	0.767	0.917		-	-	-	-	-		-	-	-	-	-
glcm-20	0.700	0.862	0.805	0.730	0.898		-	-	-	-	-		-	-	-	-	-
	2009
3-10	0.756	0.879	0.869	0.837	0.903		0.763	0.881	0.875	0.826	0.929		-	-	-	-	-
3-15	0.735	0.868	0.859	0.820	0.901		0.802	0.902	0.892	0.876	0.908		-	-	-	-	-
3-20	0.744	0.872	0.865	0.816	0.919		0.781	0.891	0.883	0.843	0.928		-	-	-	-	-
4-10	0.773	0.887	0.877	0.848	0.908		0.770	0.886	0.873	0.867	0.879		-	-	-	-	-
4-15	0.757	0.879	0.870	0.833	0.909		0.786	0.894	0.883	0.864	0.904		-	-	-	-	-
4-20	0.739	0.869	0.862	0.815	0.914		0.777	0.889	0.880	0.850	0.912		-	-	-	-	-
5-10	0.764	0.883	0.872	0.848	0.897		0.780	0.890	0.881	0.853	0.911		-	-	-	-	-
5-15	0.757	0.879	0.870	0.838	0.904		0.773	0.887	0.878	0.849	0.909		-	-	-	-	-
5-20	0.728	0.864	0.854	0.822	0.887		0.756	0.879	0.869	0.836	0.905		-	-	-	-	-
6-10	0.767	0.884	0.875	0.842	0.910		0.782	0.892	0.881	0.862	0.902		-	-	-	-	-
6-15	0.752	0.877	0.866	0.841	0.892		0.772	0.887	0.876	0.858	0.894		-	-	-	-	-
6-20	0.747	0.874	0.863	0.838	0.889		0.762	0.882	0.871	0.847	0.896		-	-	-	-	-
glcm-10	0.791	0.896	0.888	0.852	0.927		0.816	0.909	0.899	0.890	0.908		-	-	-	-	-
glcm-15	0.776	0.889	0.878	0.857	0.901		0.796	0.899	0.888	0.875	0.901		-	-	-	-	-
glcm-20	0.768	0.885	0.873	0.856	0.892		0.793	0.897	0.886	0.875	0.898		-	-	-	-	-
	2012
3-10	0.578	0.785	0.784	0.692	0.903		0.569	0.779	0.784	0.676	0.934		-	-	-	-	-
3-15	0.672	0.836	0.825	0.763	0.897		0.688	0.844	0.832	0.779	0.892		-	-	-	-	-
3-20	0.653	0.828	0.808	0.776	0.843		0.698	0.849	0.837	0.782	0.902		-	-	-	-	-
4-10	0.556	0.774	0.774	0.680	0.897		0.550	0.768	0.777	0.663	0.939		-	-	-	-	-
4-15	0.650	0.825	0.812	0.756	0.876		0.585	0.790	0.785	0.702	0.891		-	-	-	-	-
4-20	0.581	0.794	0.762	0.758	0.766		0.656	0.830	0.809	0.785	0.835		-	-	-	-	-
5-10	0.537	0.762	0.770	0.660	0.923		0.581	0.786	0.789	0.684	0.933		-	-	-	-	-
5-15	0.661	0.831	0.817	0.764	0.878		0.602	0.799	0.793	0.713	0.893		-	-	-	-	-
5-20	0.636	0.819	0.800	0.765	0.838		0.658	0.832	0.808	0.793	0.824		-	-	-	-	-
6-10	0.598	0.795	0.795	0.698	0.923		0.608	0.801	0.798	0.710	0.910		-	-	-	-	-
6-15	0.644	0.824	0.803	0.774	0.835		0.481	0.732	0.745	0.630	0.912		-	-	-	-	-
6-20	0.684	0.844	0.824	0.801	0.848		0.684	0.845	0.820	0.821	0.818		-	-	-	-	-
glcm-10	0.701	0.849	0.843	0.762	0.944		0.613	0.802	0.805	0.698	0.950		-	-	-	-	-
glcm-15	0.675	0.836	0.831	0.745	0.939		0.643	0.819	0.816	0.726	0.932		-	-	-	-	-
glcm-20	0.659	0.830	0.816	0.764	0.876		0.597	0.802	0.772	0.765	0.780		-	-	-	-	-
	2015
3-10	0.414	0.711	0.763	0.661	0.902		0.549	0.777	0.813	0.716	0.941		0.704	0.853	0.870	0.800	0.955
3-15	0.372	0.691	0.758	0.372	0.936		0.452	0.730	0.783	0.452	0.943		0.698	0.850	0.870	0.698	0.970
3-20	0.321	0.667	0.748	0.614	0.957		0.467	0.738	0.790	0.672	0.958		0.691	0.847	0.868	0.782	0.974
4-10	0.422	0.715	0.764	0.667	0.894		0.498	0.752	0.794	0.695	0.925		0.678	0.840	0.860	0.786	0.950
4-15	0.408	0.708	0.763	0.656	0.913		0.526	0.766	0.805	0.705	0.938		0.648	0.826	0.852	0.758	0.973
4-20	0.431	0.719	0.774	0.662	0.931		0.483	0.745	0.795	0.679	0.959		0.656	0.830	0.856	0.761	0.977
5-10	0.388	0.698	0.749	0.655	0.873		0.469	0.737	0.780	0.687	0.902		0.668	0.835	0.857	0.775	0.959
5-15	0.433	0.721	0.773	0.666	0.921		0.518	0.762	0.801	0.704	0.930		0.647	0.825	0.851	0.761	0.965
5-20	0.434	0.721	0.776	0.662	0.937		0.494	0.751	0.799	0.683	0.963		0.623	0.813	0.844	0.743	0.976
6-10	0.416	0.713	0.772	0.654	0.942		0.491	0.749	0.797	0.683	0.958		0.669	0.836	0.860	0.768	0.977
6-15	0.402	0.706	0.767	0.648	0.940		0.489	0.748	0.796	0.684	0.953		0.634	0.819	0.847	0.752	0.970
6-20	0.377	0.694	0.763	0.635	0.956		0.485	0.746	0.797	0.679	0.963		0.642	0.823	0.848	0.759	0.962
glcm-10	0.398	0.704	0.770	0.643	0.961		0.593	0.798	0.827	0.741	0.935		0.691	0.847	0.867	0.785	0.967
glcm-15	0.246	0.631	0.730	0.587	0.967		0.425	0.717	0.778	0.654	0.962		0.709	0.856	0.874	0.794	0.974
glcm-20	0.264	0.640	0.736	0.592	0.970		0.397	0.704	0.770	0.643	0.959		0.591	0.798	0.833	0.727	0.976

and in the diagrams in Figure 4.

Generally, it was found that the best results for most archival materials were obtained using a local texture analysis radius of 10 or 15 pixels. A larger radius of the area of analysis (20 pixels) gave better results only in 2012 (early spring). Visualizations of the results of the texture analysis operation at various analysis parameters (before closing) for a fragment of the research area are presented in the figures included in Supplementary 1 (Figures S1–S14).

In the case of black and white aerial photographs (P) from 1971, the highest accuracy was obtained using the GLCM method and a radius of texture analysis of 10 pixels (Table 4, Figure 4). Slightly lower values of individual accuracy parameters were obtained for granulometric analysis p-1971-4-10, p-1971-3-10, and p-1971-5-10.

For 1996, the best results were obtained based on granulometric analysis: c-1996-3-15, c-1996-3-10, p-1996-3-15, and p-1996-3-10. The GLCM method for these photographs is characterized by the lowest values of all analyzed accuracy parameters (Table 4, Figure 4).

In the case of 2003, relatively small differences were observed in the evaluation of the accuracy of individual variants of granulometric analysis and GLCM (Kappa coefficient from 0.70 to 0.76) (Table 4). The variants with the highest accuracy were: p-2003-3-10, p-2003-4-15, p-2003-5-15, and p-2003-6-15.

In turn, as far as 2009 is concerned, the best results were obtained using the GLCM method, and this applies to the analysis of both the simulated image P and the original RGB images (Table 4, Figure 4). The highest accuracy was obtained for the following variants: c-2009-glcm-10, c-2009-3-15, and p-glcm-2009-10.

The GLCM method proved to be the most effective in relation to the imagery from 2012—the variant p-2012-glcm-10 was characterized by the highest accuracy (OA, kappa, F1, UA) (Table 4, Figure 4). In the case of granulometric analysis, the best results were obtained by analyzing RGB images in the variant c-2012-3-20.

Imagery from 2015 enabled to carry out the widest list of experiments—using the simulated P picture, as well as RGB and CIR images. Definitely the highest accuracy was obtained on the basis of CIR image analysis (Table 4, Figure 4). The highest values of individual accuracy indicators were characterized by the following variants: cir-2015-glcm-15, cir-2015-3-10, cir-2015-3-15, and c-2015-3-20.

Comparing the results obtained on the basis of RGB images and of the P images simulated on their basis (Figure 4) in most cases, higher accuracy was obtained for RGB images.

6. Discussion

The tested variants differed in many features: Type of spectral data (P, RGB, CIR), type of textural data (method and selected processing parameters), but also in the quality of the source spectral data related to the type of images (analog, digital), original image scale, original spatial resolution, and general radiometric quality of photos. However, the analysis of the results allows one to see some patterns.

Firstly, the comparison between variants differing only with the type of spectral data used in the classification (P/RGB/CIR; images from 1996, 2009, 2012, 2015) unambiguously indicate a strong relationship between spectral resolution or the type of spectral bands used, and the accuracy of classification, also using textural analysis. Variants using RGB images were more accurate compared to P in almost all cases (Table 4, Figure 4). The only case in which this relationship is not completely unambiguous is the classification of the photo from 2012; however, this is an early spring photo, and at the time of exposure, the deciduous plants were free of (visible) leaves (see Figure 3), which caused the spectral differences between them and the other classes of land cover to be much smaller than in the case of photos from other periods (with visible leaves). This affected both the texture of wooded and bushy areas, as well as the areas of meadows that dominate in this area.

In addition, the case of the imagery from 2015 confirmed the thesis about the importance of the near-infrared range (Figure 5, Figure 6). The accuracy of classification of selected variants of this imagery was generally low (measured by the Cohen’s kappa coefficient, the F1-score achieves quite high values); however, variants using CIR images were significantly better than the other variants: The kappa coefficient values were at 0.1–0.2 higher than for RGB variants and 0.2–0.3 higher than for P variants (Table 4). Also, PA and UA are the highest in the case of CIR image analysis (Figure 5). The difference in accuracy (measured by the Cohen’s kappa coefficient and the F1-score) between the P and RGB variants is lower than between RGB and CIR, and this is the case for each analyzed variant of textual analysis (Figure 5). It shows the important role of near-infrared range in the analysis of vegetation. This spectral band, as more sensitive to the diversity of plant cover, allows better detection of trees and shrubs. In the case of P and RGB photographs, there is an overestimation of the ranges of trees and shrubs for herbaceous vegetation (Figure 6), with the largest overestimation for the P imagery.

The above conclusions correspond with the results obtained in previous studies, showing the importance of appropriate spectral resolution in the classification using both spectral data and texture analysis results [68,69], and the importance of image quality for the results of texture analysis itself [70].

Another basically unambiguous conclusion is the usefulness of the proposed additional processing that is closing by reconstruction, applied to individual masks of wooded areas obtained as a result of classification. Closing is an extensive operation, which means that in its result, the mask can only increase (or remain unchanged). At the same time, the use of the operation by reconstruction does not lead to a change in the shape of the original mask, but only to filling of relatively small holes. According to a detailed accuracy analysis, the use of the closing by reconstruction operation led to a decrease in UA, but at the same time, to a more significant increase in PA and, as a result, to an increase in the overall accuracy of classification. This is presented in Figure 7. Unfortunately, in the case of small trees and shrubs growing in a large dispersion, this may lead to overestimation of the area occupied by them, which is also seen in the analysis of texture for 2012 presented for two variants (granulometry and GLCM) in Figure 8.

The comparison of the effectiveness of two methods: GLCM and granulometric analysis do not give a clear indication of the better one. In two cases (1971 and 2009), better results were obtained by using the GLCM method, although the differences between the two methods were relatively small (the difference between kappa values between the best variants was around 0.05 for 1971 and 0.01 for 2009, in favor of GLCM). In one case (2015), the accuracy obtained for the best variants measured by the kappa coefficient and other general indicators were basically identical; the difference in favor of one of the GLCM variants was 0.005, which can be considered as negligible. In the remaining three cases, greater accuracy was observed for variants based on granulometric analysis. In two cases (2003 and 2012), the differences in accuracy were very small (differences between the kappa coefficient values for the best variants amounted to approximately 0.01, in favor of granulometric analysis). However, in the case of 1996, these were very large differences (between 0.2 and 0.5) in favor of granulometric analysis, which is mainly due to the very low accuracy of GLCM based variants. This is a case worth a separate commentary.

Analyzing Figure 4 in the context of the dates of imagery acquisition, it can be seen that the best results for photos from the period of dense foliage (May–September, with the exception of data from 1971 with lower radiometric quality) gave the use of granulometric analysis using six granulometric maps and a radius of 10–15 pixels (variants marked as x-xxxx-3-10 and x-xxxx-3-15). However, in the case of photos from the spring period (March–April), with a less dense foliage, the highest accuracy was obtained by using the GLCM method with a radius of 10–15 pixels.

The masks obtained on the basis of the GLCM analysis or granulometric analysis are presented in Figure 9. The low accuracy of variants based on the GLCM analysis results mainly from very low user accuracy. This is due to the high texture of agricultural areas (resulting from agrotechnical activities, including mowing meadows) and the areas on which the trees were cleared. At the same time, granulometric analysis results in a much smaller error of commission (Figure 9, Table 4). Probably, this is due to its multi-scale: Areas characterized by high but different textures are marked differently on different granulometric maps. This increases the mutual distinction of these areas. In addition, the accuracy of the GLCM-based variants for 1996 clearly decreases as the area of the local texture analysis increases. For an area with a radius of 10 pixels, the kappa value is 0.597, for a radius of 15 pixels is 0.572, while for an area with a radius of 20 pixels is only 0.305. As can be seen in Table 4, the decrease in the overall accuracy is mainly caused by a decrease in PA, and thus an increase in the error of omission. Analysis of the images of individual GLCM indices shows that part of the tree-covered areas obtains identical values as unripen areas. The area of these types of surfaces increases with the radius of the analysis area. This can be seen in Figure 9: The boundary of the wooded area detected on the basis of the selected GLCM indices “moves back” inside the area, along with the increase of the analysis area. In the case of granulometric analysis, this effect is not so pronounced, although it also occurs. Here too, the importance of multiscale granulometric analysis can be seen. Different values in selected granulometric maps indicate the “activation” in response to a texture with grains of different sizes. In GLCM indicator images, areas with a larger grain texture simply resemble areas with no clear texture. Hence the observed decrease in accuracy for the classification based on the GLCM analysis.

Interestingly, in the case of other images from different dates, this effect is not so important. It seems that this can be explained by the relatively low (in some areas) clearness of the texture in the 1996 picture compared to, for example, the picture from 2009 (compare Figure 3), which in turn may result from the small scale of the original image from 1996, 1:26 000.

Based on the above, it can be concluded that granulometric analysis is more stable and, in general, slightly better as a texture descriptor for image classification. There is no extensive literature on this subject, but existing publications also show greater efficacy of granulometric analysis in this area [36,71].

The analysis of the influence of the size of a texture analysis area on the accuracy of the classification indicates a general (though not unequivocal) tendency for the classification quality to decrease along with the increase of the size of the analysis environment (figures for all variants may be found in the Supplementary Materials). In most cases, the worst results were obtained for the analysis performed for a 20-pixel environment. The only exception is the classification of images from 2012; these, however, differ significantly from the others when it comes to the image of wooded and bushy areas (early spring photo and lack of leaves on deciduous trees and shrubs). However, when analyzing areas covered with coniferous trees, a similar tendency can be noticed as described above—a decrease in accuracy when the area of analysis increases. One can also observe a certain consistency in terms of individual variants within one source image—the relationships between variants differing only in the size of area of a texture analysis are similar for different texture analysis methods and for a different number of granulometric maps used (in the case of the granulometric analysis).

A comparison of the results of classification variants based on a granulometric analysis (differing in the number of granulometric maps used) allows to draw unambiguous conclusions. For this type of texture and with this pixel size (0.5 m), the use of granulometric maps based on a structuring element with a radius of more than 3 pixels (a size of approximately 3.5 m) does not improve the classification efficacy, and in most cases on the contrary: The use of granulometric maps showing the presence of textures with larger grains (larger than 3.5 m) can lead to a decrease in efficacy. This is due to the fact that this type of information is no longer related to the occurrence of wooded or shrubbed areas.

Analyzing the course of the border of wooded and bushy areas for particular dates and variants of texture analysis, the impact of the date of image acquisition is clearly visible. The problem is the shadow, which in the pictures acquired in early spring (2012) is the longest and increased significantly the detected surface of wooded areas (Figure 10). On the other hand, the shadow increases the texture of the image of deciduous trees—leafless at this time of the year. Due to this effect, they were correctly classified. In other dates, this effect is weaker—the borders were best reproduced in the photographs from August (2015) and May (1996 and 2003), i.e., the period when the sun is highest above the horizon, and the shadows are short.

The analysis of the results also shows that regardless of the archival data and the texture analysis method used, individual trees and shrubs are difficult to detect—only larger trees and shrubs are possible to detect with these methods.

The importance of impact of the scale/GSD of aerial photographs on the effectiveness of determining the range of trees and shrubs has not been observed. The most important were the impact of the time period and radiometric quality of the source materials. Some of the errors in the case of photos from 1971, including film damage (scratches were seen in the images), resulted from the poor quality of these sources.

Summing up, in general, the highest accuracy in determining the range of trees and shrubs for photos recorded in the summer was obtained using granulometric analysis with six granulometric maps and a analysis radius of 10–15 pixels. For photos from the spring period, the GLCM method gave better results, also with a 10–15 pixel radius of analysis.

When assessing the tested methods from the point of view of their possible application, it should be noted that they are easy to implement and automate the classification process (also in open-source software). Importantly, compared to the dense image matching (DIM) technique (one of the most effective methods of analyzing aerial photographs in the context of research), they enable the determination of the range of trees and shrubs based on images obtained during the period when the trees and shrubs are not fully covered by leafs (early spring). DIM, although a method characterized by high accuracy in determining object boundaries and the ability to analyze in three dimensions, does not allow the determination of the range of trees and shrubs in the leafless period [20]. During this time, it is only possible to analyze the range of coniferous trees and shrubs, which may be of interest, e.g., in the context of conducting research on the impact of climate change on the range of occurrence of individual species of trees and shrubs [72].

7. Conclusions

The aim of the research was to assess the applicability of texture analysis for the analysis of the dynamics of succession of trees and shrubs. The research was carried out using aerial imagery (characterized by various technical parameters) acquired in six different years (1971, 1996, 2003, 2009, 2012, and 2015) in various phenological periods.

A total number of 330 classification variants, which differ in source data (date of image, spectral resolution) and texture data (method and parameters of texture analysis) were tested. The analysis of the obtained results showed great effectiveness in determining the wooded and shrubbed areas using texture analysis.

In the majority of cases tested, the results obtained indicate a similar efficacy of both methods of texture analysis (in the analyzed range), except for one case (1996), in which the efficacy of granulometric analysis was significantly better than the efficacy of GLCM analysis, due to the high texture in areas of arable land and meadows and to the relatively small scale of the photos. The advantage of the granulometric analysis resulted mainly from multi-scality, thanks to which the wooded/shrubbed areas were distinguished better from agricultural areas with a relatively high texture. Therefore, it appears from all the analyzed examples that granulometric analysis is much more stable in terms of the analyzed application. This confirms the results of previous studies which compared, among others, these two above-mentioned methods [36,71].

The accuracy of determining the wooded and shrubbed areas is also influenced by the date of image acquisition—in the case of early spring images, the observed long shadows, typical for this period, caused a significant extension of the area of tree/shrub masks. However, the role of the shadows was also positive, because thanks to them, the texture of the areas with deciduous trees, which were leafless at this time of year, was strengthened. It is therefore possible to designate tree and shrubbed areas based on images from leafless periods using texture analysis, which is much more difficult when using dense image matching methods [20].

In addition, in summer time, wooded and shrubbed areas may be overestimated due to the presence of herbaceous plants, which can also be characterized by high textures.

The operation of closing by reconstruction proved to be an effective method of filling gaps in wooded areas. However, in the case of large scattering of trees and shrubs, it leads to overestimation of the range of such areas. Therefore, it is necessary to consider the legitimacy of its use in areas where the succession of trees and shrubs can be thus characterized.

Summing up the conducted analyses, it can be concluded that texture analysis is an effective method for determining the range of wooded and shrubbed areas, in particular those with considerable density of trees and shrubs. The granulometric analysis showed generally greater suitability for this purpose than GLCM. In the case of detection of individual trees and shrubs, the effectiveness of such analyses is smaller—large trees and shrubs are generally correctly detected, while small-sized objects are not. This means that these methods can be an effective tool for monitoring only the later stages of succession of trees and shrubs.