**4. Materials and Methods**

The main stages of the development methodology included:

	- a. Selection of methods and parameters of texture analysis;
	- b. Creation of spectro-textural data sets;
	- c. Classification using support vector machine (SVM); and
	- d. Additional processing.

The methodology scheme is presented in Figure 2.

**Figure 2.** The methodology scheme.

### *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, and in Figure 3, selected parts of photos are visualized.


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

**Figure 3.** Parts of the aerial photos used for analysis.

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 (*Ps*) were calculated by averaging the values of red (*R*) and green (*G*) bands:

$$P\_s = \frac{R+G}{2},\tag{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:


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:


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.


**Table 2.** Tested variants of data sets.

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


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

#### 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 m2; 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 and in the diagrams in Figure 4.


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


**Table 4.** *Cont.*

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.

**Figure 5.** Comparison of accuracy indicators (Cohen's kappa, F1-score, Producer's Accuracy—PA, User's Accuracy—UA) for the classification based on individual granulometric analysis variants obtained for 2015.

**Figure 6.** Comparison of tree/shrub masks obtained on the basis of images from 2015: P, RGB, and CIR in two variants of texture analysis—granulometry and GLCM.

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.

**Figure 7.** Comparison of tree and shrub masks for 2009, variant c-2009-3-15 before closing and after closing, with a reference mask imposed.

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.

**Figure 8.** Comparison of tree and shrub masks for 2012, variants c-2012-3-10 and c-2012-glcm-10 before closing and after closing, with a reference mask imposed.

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

**Figure 9.** Comparison of tree and shrub masks for 1996, variants c-1996-3-10, c-1996-3-15, c-1996-3-20, c-1996-glcm-10, c-1996- glcm-15, and c-1996-glcm-20, with a reference mask imposed.

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.

**Figure 10.** The effect of the date of image acquisition on the accuracy of determining the boundaries of wooded and shrubbed areas.

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.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2220-9964/8/10/450/s1: Figure S1: Comparison of tree and shrub masks for 1971 obtained using P image, granulometric analysis variants, with a reference mask imposed; Figure S2: Comparison of tree and shrub masks for 1996 obtained using P image, granulometric analysis variants, with a reference mask imposed; Figure S3: Comparison of tree and shrub masks for 1996 obtained using RGB image, granulometric analysis variants, with a reference mask imposed; Figure S4: Comparison of tree and shrub masks for 2003 obtained using P image, granulometric analysis variants, with a reference mask imposed; Figure S5: Comparison of tree and shrub masks for 2009 obtained using P image, granulometric analysis variants, with a reference mask imposed; Figure S6: Comparison of tree and shrub masks for 2009 obtained using RGB image, granulometric analysis variants, with a reference mask imposed; Figure S7: Comparison of tree and shrub masks for 2012 obtained using P image, granulometric analysis variants, with a reference mask imposed; Figure S8: Comparison of tree and shrub masks for 2012 obtained using RGB image, granulometric analysis variants, with a reference mask imposed; Figure S9: Comparison of tree and shrub masks for 2015 obtained using P image, granulometric analysis variants, with a reference mask imposed; Figure S10: Comparison of tree and shrub masks for 2015 obtained using RGB image, granulometric analysis variants, with a reference mask imposed; Figure S11: Comparison of tree and shrub masks for 2015 obtained using CIR image, granulometric analysis variants, with a reference mask imposed; Figure S12: Comparison of tree and shrub masks for 1971, 1996, and 2003, GLCM variants, with a reference mask imposed; Figure S13: Comparison of tree and shrub masks for 2009 and 2012, GLCM variants, with a reference mask imposed; Figure S14: Comparison of tree and shrub masks for 2015, GLCM variants, with a reference mask imposed.

**Author Contributions:** The aim of the study and the methodology of research: P.K. and K.O.-S.; texture analysis and classification: P.K. and K.L.; reference data preparation: A.P. and K.O.-S.; accuracy assessment: K.O.-S.; results analysis, preparation of tables and charts, synthesis of the study data: K.O.-S. and P.K.; writing—original draft preparation: P.K. and K.O.-S.; writing—review and editing: K.O.-S. and P.K.; visualization preparation: K.O.-S.; WP5 leading and supervision: K.O.-S.

**Funding:** This study was co-financed by the Polish National Centre for Research and Development (NCBiR), project No. DZP/BIOSTRATEG-II/390/2015: The innovative approach supporting monitoring of non-forest Natura 2000 habitats, using remote sensing methods (HabitARS). The Consortium Leader is MGGP Aero. The project partners include: University of Lodz, University of Warsaw, Warsaw University of Life Sciences, Institute of Technology and Life Sciences, University of Silesia in Katowice, and Warsaw University of Technology. An article processing charge was financed from the statutory subsidy of the Faculty of Geodesy and Cartography, Warsaw University of Technology.

**Acknowledgments:** We would like to express our sincere gratitude to Bo ˙zena Michna and Agnieszka Jakubowska for administrative support during HabitARS project.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analysis, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

## **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*ISPRS International Journal of Geo-Information* Editorial Office E-mail: ijgi@mdpi.com www.mdpi.com/journal/ijgi

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18