**Przemysław Kupidura \*, Katarzyna Osi ´nska-Skotak, Katarzyna Lesisz and Anna Podkowa**

Department of Photogrammetry, Remote Sensing and Spatial Information Systems, Faculty of Geodesy and Cartography, Warsaw University of Technology, 1, 00-661 Warsaw, Poland; katarzyna.osinska-skotak@pw.edu.pl (K.O.-S.); lesisz.katarzyna@gmail.com (K.L.);

anna.podkowa.dokt@pw.edu.pl (A.P.)

**\*** Correspondence: przemyslaw.kupidura@pw.edu.pl; Tel.: +48-22-234-5389

Received: 20 August 2019; Accepted: 10 October 2019; Published: 12 October 2019

**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.

**Keywords:** secondary succession monitoring; Natura 2000 threats; tree detection; archival photographs; spectro-textural classification; granulometric analysis; GLCM

#### **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–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–13]. Sometimes, stereodigitalization is also used [14–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–23], fractal analysis [24], discrete wavelet transformation [25], Laplace filters [26–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–24,33–36,38–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–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]).

$$\text{Energy} = \sum\_{\mathbf{i}, \mathbf{j}} \mathbf{g}(\mathbf{i}, \mathbf{j})^2 \tag{1}$$

Energy is a description of the uniformity of texture.

$$Entropy = \sum\_{i,j} g(i,j) \log\_2 g(i,j), \text{ or } 0 \text{ if } g(i,j) = 0,\tag{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 = \sum\_{i,j} \frac{(i-\mu)(j-\mu)g(i,j)}{\sigma^2} \tag{3}$$

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

$$\text{Inverse Difference Moment} = \sum\_{i,j} \frac{1}{1 + \left(i - j\right)^2} g(i, j), \tag{4}$$

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

$$Inertia = \sum\_{i,j} (i,j)^2 g(i,j)\_{\prime} \tag{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.

$$\text{Cluster Shade} = \sum\_{i,j} ((i - \mu) + (j - \mu))^3 \lg(i, j), \tag{6}$$

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

$$\text{Cluster Premience} = \sum\_{i,j} \left( (i - \mu) + (j - \mu) \right)^4 \text{g}(i, j), \tag{7}$$

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

$$\text{Haralick's Correlation} = \frac{\sum\_{i,j}(i,j)g(i,j) - \mu\_t^2}{\sigma\_t^2},\tag{8}$$

where (*i*, *j*) is the matrix cell index, *g*(*i*, *j*) is the frequency value of the pair having index (*i*, *j*), μ = *<sup>i</sup>*,*<sup>j</sup> i*∗ *g*(*i*, *j*) = *<sup>i</sup>*,*<sup>j</sup> j*∗ *g*(*i*, *j*) (due to matrix symmetry) and means weighted pixel average, σ = *<sup>i</sup>*,*j*(*i* − μ) <sup>2</sup> <sup>∗</sup> *<sup>g</sup>*(*i*, *<sup>j</sup>*) <sup>=</sup> *<sup>i</sup>*,*j*(*j* − μ) <sup>2</sup> <sup>∗</sup> *<sup>g</sup>*(*i*, *<sup>j</sup>*)(due to matrix symmetry) and means weighted pixel variance, and μ*<sup>t</sup>* and σ*<sup>t</sup>* 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].

**Figure 1.** Location of the study area (blue rectangle) and overview (prepared based on: http: //geoserwis.gdos.gov.pl/mapy/, photo: E. Sierka). The photo shows a fragment of valuable non-forest habitats, where secondary succession is observed.

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 km2) 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.
