Assessment of Textural Differentiations in Forest Resources in Romania Using Fractal Analysis

Abstract

Deforestation and forest degradation have several negative effects on the environment including a loss of species habitats, disturbance of the water cycle and reduced ability to retain CO₂, with consequences for global warming. We investigated the evolution of forest resources from development regions in Romania affected by both deforestation and reforestation using a non-Euclidean method based on fractal analysis. We calculated four fractal dimensions of forest areas: the fractal box-counting dimension of the forest areas, the fractal box-counting dimension of the dilated forest areas, the fractal dilation dimension and the box-counting dimension of the border of the dilated forest areas. Fractal analysis revealed morpho-structural and textural differentiations of forested, deforested and reforested areas in development regions with dominant mountain relief and high hills (more forested and compact organization) in comparison to the development regions dominated by plains or low hills (less forested, more fragmented with small and isolated clusters). Our analysis used the fractal analysis that has the advantage of analyzing the entire image, rather than studying local information, thereby enabling quantification of the uniformity, fragmentation, heterogeneity and homogeneity of forests.

1. Introduction

Forests have been managed intensively due to increased demand for wood and cropland [1,2,3]. Sustainable forest management is imperative when facing important environmental and socio-economic changes in forest resources, as pressure on the forests’ resources has an adverse effect on several components of ecosystems (water pollution, reducing the biodiversity, climate change, soil erosion) [4]. Specifically, long-term changes in forest environments have direct effects on climate change, carbon cycling, and carbon–nitrogen interactions, thus underlining the need for novel management strategies to maintain and preserve the global carbon storage and ecosystem roles of the forest [5,6]. The European continent hosts a variety of forest biomes, from the temperate to boreal regions, generally characterized by considerable anthropogenic influence causing environmental changes [7,8,9,10].

Romania incorporates large continuous virgin forest ecosystems in the temperate zone of Europe [11,12]. In the recent decades, the deforestation rate in Romania has rapidly intensified in the Carpathians and Subcarpathians through a surge in the demand from the wood industry, timber export growth and manufacturing of wood [3,13]. The natural forest cover in Romania is under threat due to the high population growth rate [14]. Urbanization processes (new buildings for residential, commercial, academic and business purposes) and the conversion of rural areas into urban developed areas have caused substantial changes in land use and land cover [15].

Romania is a competitive player on the European wood market. It produces a quarter of European hardwood lumber output due to it having the lowest selling prices on the continent [16]. According to the National Institute of Statistics [17], the amount of roundwood produced in 2002 was 4.3 million m³ from coniferous forests, 4 million m³ from broadleaved forests, 1.1 million m³ from hardwood forests and 1.1 million m³ from softwood forests. A considerable increase in production has happened since the turn of the new millennium with 2015 values of 6 million m³ coniferous forest, 6.6 million m³ broadleaved forest, 5 million m³ hardwood forest, and 1.1 million m³ softwood forest. The Carpathian Mountains generally have a total growing stock of about 250–300 m³/ha for coniferous species with moderately lower the stock rates in the Western Carpathian Mountains (about 200–250 m³/ha) [18]. The utilization rate of forest decreased to 48% in 2005 expressed as a percentage of the annual increment, a decrease from 60% in 1990 [19].

Based on the forest management literature in Romania, sustainable forestry should consider:

• The preventive role of forests through the stability of slopes mitigating landslides, floods and flash floods, and catastrophic inundations of river watersheds [20,21,22,23].
• The ecological function of forests (providing ecosystem services such as biodiversity) [24].
• The socio-economic trends governed by political changes after the post-socialist period [25,26] and forest restitutions [27,28] to former private owners and legal entities [29].

The assessment of forest cover change is a first necessary step to identify critical areas of land use changes [30]. Methods of forest change mapping are mainly based on an assessment from remotely sensed imagery based on changes in spectral characteristics though time [30,31] or from indices designed for mapping forest changes such as the disturbance index [32,33,34]. This index was applied to map forest cover changes in the protected areas of the Eastern Carpathians Mountains in Romania (Maramureş Mountains Nature Park, Rodna Mountains National Park, Călimani National Park) [28].

However, a spatially explicit identification and quantification of forest fragmentation, together with the spatial patterns and dynamics quantification during deforestation processes are needed. The powerful quantitative approach for this was presented by [35] based on the fractal analysis (fractal geometry). The fractal analysis provides a complementary approach to traditional per-pixel based assessment of deforestation and forest disturbance, by including a contextual dimension into the analysis. Fractal analysis in the form of the box-counting, pyramid dimension, Minkowski dimension, Tug-of-War dimension and Fractal Fragmentation Index (FFI) were recently used in the analysis of forests [36,37,38].

In this study, we used a different approach to fractal analysis in relation to textural differentiations of forest resources by analysis of the fractal object dilation effects on fragmentation and heterogeneity of the forest spatial evolution. We applied this approach to forest resources from the eight development regions (DRs) in Romania to obtain information about the spatial effects and distribution of deforestation and reforestation. Through fractal analysis, we provide a quantification of the degree of fragmentation of forested areas. Ultimately, this research contributes by assisting institutional interventions to alleviate the imbalances that are the consequence of deforestation in Romania.

2. Materials and Methods

2.1. Study Area

This study includes eight administrative regions (development regions; DRs) in Romania, namely, North-East, South-East, South Muntenia, South-West Oltenia, West, North-West, Center and Bucharest-Ilfov DRs (Figure 1).

The development regions of Romania were established in 1998 to facilitate the effective regional coordination in Romania following the European Union (EU) accession. According to the law 315/2004, development regions are not administrative-territorial units and are therefore not legal entities of the counties, thus posing challenges to forest management [39]. Forest cover in Romania of 63,990 km² includes beech (Fagus sylvatica), fir (Abies alba), spruce (Picea abies), oak (Quercus robur), ash (Fraxinus excelsior) and pine (Pinus) [40]. The largest forest area of 51.9% is in the Carpathian Mountains (formed by coniferous and deciduous species) followed by 37.2% in the hills (the Subcarpathians, the Getic Plateau, the Transylvanian Plateau, the Western Hills) and 10.9% in the plains (the West Plain and the Romanian Plain). Forest exploitations in Romania is in continuous growth in the Carpathian Mountains and in the hilly areas [20]. The most deforested areas are in the Eastern Carpathians (counties Suceava (SV), Maramureș (MM), Harghita (HR), Covasna (CV), Vrancea (VN)) and in the Southern Carpathians (counties Argeș (AG), Vâlcea (VL), Hunedoara (HD)).

The development regions with the most intense deforestation are the regions with the largest forest areas, characterized by a mountain relief and high hills: North-East Region (599.31 km²), Center Region (581.19 km²) and North-West Region (505.04 km²). The regions with the lowest forest cover, dominated by low hills and plains, are generally less deforested (the West Region (194.57 km²), South Region (139.71 km²), South-East Region (98.39 km²) and Bucharest-Ilfov Region (3.16 km²) [40]. The most intense deforestation has occurred in North-East and North-West Regions (deforestation rates of 4.8% and 3.9% respectively), followed by the Center Region (3.9%), South Region (2.5%), South-West Region (1.5%), West Region (1.4%), South-East Region (1.4%) and Bucharest-Ilfov Region (1.2%) [36].

The analysis was conducted at the level of development regions because of the pronounced differences in the way forest resources have been managed between regions which in turn reflected the levels of deforestation and reforestation.

2.2. Image Preprocessing

Digital images obtained from the GFC (Global Forest Change) 2000–2012 database, the Department of Geographical Sciences at the Maryland University were used as input for the fractal analysis. This database is the result of analyzing Landsat 7 ETM+ (30 m spatial resolution) data of forest areas during the period 2000–2012 [41], representing global information on forest disturbance with the highest spatial resolution currently available. The Landsat spatial resolution provides sufficient spatial detail for a good fractal analysis and we used a time-series of satellite images representing the annual time steps from 2000 to 2012 to quantify and analyze areas of deforestation and forest reforestation.

The initial projection of the Landsat GFC data was converted from WGS 1984 World Mercator (EPSG 3395) to the Stereographic 70 coordinate system (EPSG 31700). Specific information about Landsat imagery used for the analysis mask is presented in

Table 1. Landsat 7 ETM+ image paths, rows

Satellite images	Resolution	Longitude	Latitude	Paths	Rows
LANDSAT 7 ETM+	30 m	50°0’ N	20°18’ E	188	25
LANDSAT 7 ETM+	30 m	50°0’ N	30°18’ E	181	25
LANDSAT 7 ETM+	30 m	39°60’ N	20°18’ E	185	32
LANDSAT 7 ETM+	30 m	39°60’ N	30°18’ E	179	32

. The satellite images were stored in uncompressed color TIFF format by using ArcGIS and the raster images were transformed to vector (Raster to point) with the resolution of each pixel being equivalent to the dimension of a point. To extract surfaces for each territorial-administrative unit (UAT), a spatial join was conducted, allowing for points to be grouped for each territorial-administrative unit. The calculation of forest change areas for each UAT and for each year from 2000 to 2012 was repeated for the three types of Landsat 7 ETM+ assessments available in the GFC database (forest cover loss, forest cover gain and year of gross loss).

The Landsat 7 ETM+ imagery was segmented using the Color Deconvolution vector H & E DAB plugin for the ImageJ software version 1.49t [42] and subsequently converted to binary images, using a threshold of 1–255. The binarized images only contained black pixels with 0 value (background pixels, non-fractal pixels, meaning non-forest areas) and white pixels with 255 values (foreground pixels, fractal pixels, meaning forested, deforested and reforested areas).

2.3. Fractal Analysis

We conducted a fractal analysis of the digital binarized images of the forest area evolution affected by deforestation and reforestation for eight DRs. Fractal analysis was conducted because of the ability to provide additional textural information on the effects incurred by deforestation which provides complementary information on forest disturbance as compared to traditional per-pixel forest disturbance analysis. Fractalyse 2.4. software (Research Centre ThéMA, University of Franche-Comté, Besançon, France) was employed for fractal analysis [43,44].

The box-counting dimension was used for determining the fractal dimension of the forest areas (D_Surf-dil0), the fractal dimension of the dilated forest areas (D_Surf-dil3) and the fractal dimension of the border of the dilated forest areas (D_Bord-dil3). This is the most used method to estimate fractal dimensions. The image is covered by a quadratic grid and the grid size ε is varied. For each value ε, the number of squares N(ε) containing any black pixels are counted. In Fractalyse 2.4, grids of different sizes are overlaid on a two-dimensional black and white image of forest areas and the number of black cells is counted for each grid size. The black cells represent the forest areas while the white ones indicate an absence of forest areas. A double logarithmic plot of N(ε), the number of occupied cells versus 1/ε is drawn and subsequently the slope of a linear regression gives an estimate of the box dimension.

The mathematical formula used for the calculation of the box counting fractal dimension was:

$$D_{B-C}=\underset{\epsilon \to 0}{\lim }\left(\frac{\log N\left(\epsilon \right)}{\log \frac{1}{\epsilon }}\right)$$

(1)

where D_B-C is the box-counting fractal dimension, ε is the side length of the box, and N(ε) is the number of non-overlapping and contiguous boxes of side ε required to cover the area of the fractal object [40,45]. As the zero limits cannot be applied to digital images, D_B-C was estimated by means of the equation:

$$D_{B-C}=d$$

(2)

where d is the slope of the graph log (N(ε)) against log (1/ε) [45].

The box-counting algorithm cannot offer a perfect coverage of a fractal pattern, so the number of occupied boxes is approximative. Hence, the estimation of the box-counting dimension requires a sufficient number of levels of analysis [43] but also the largest possible range of box sizes [46].

Four fractal dimensions were considered in this study: the fractal box-counting dimension of the forest areas (D_Surf-dil0), the fractal box-counting dimension of the dilated forest areas (D_Surf-dil3), the dilation fractal dimension (D_dil3) and the fractal box-counting dimension of the border of the dilated forest areas (D_Bord-dil3). We analyzed D_Surf-dil3, D_dil3 and D_Bord-dil3 to determine whether large variation in the fractal dimension was a resultof a transition from the fractal organization with a lot of free spaces (D_Surf-dil0) to another organization with less free spaces or a more compact organization (D_Surf-dil3) (Figure 2).

The first step of the morphological characterization of the forest patterns involves measuring the D_Surf-dil0 which describes how the forest mass is concentrated across scales on a given area. D_Surf-dil0 can take any value between 0 and 2. D_Surf-dil0 equal to 0 corresponds to a limited case where the pattern is made up of a single point (e.g., a single patch of forest was only 0.022 km²). D_Surf-dil0 < 1 corresponds to a pattern made up of unconnected elements, detached clusters (the patches of forests are highly dispersed and mostly isolated from one another, similarly to a Fournier Dust). D_Surf-dil0 > 1 corresponds to a mix of connected elements forming large clusters, connected elements forming small clusters, and isolated elements. When D_Surf-dil0 is close to 2, the forest pattern is quite uniformly distributed following only a single scale with a high degree of homogeneity. In this case, the forest is uniformly distributed and most elements are connected to each other. The overall shape is then mainly a single large and compact cluster. An increasing D_Surf-dil0 indicates a more dense and homogenous pattern of forest areas with increasing filling of the existing areas characterized by gaps. A decreasing D_Surf-dil0 indicates a growing fragmentation in multi-scale clusters and a more uneven distribution of forest areas. This uneven fragmentation of deforested or reforested areas can be imposed by human activity or natural conditions by themselves. The quality of the fractal dimension estimation is controlled by using the Pearson correlation coefficient R. Estimations were accepted when R exceeded the value of 0.99 [47]. If the fit between the empirical curve and the estimated curve is low, it is be concluded that either the pattern being studied is not fractal, or that it is multi-fractal [43].

The second step of the morphological characterization of the forest pattern involves a measure of the D_dil3. The principle of dilation is based on an algorithm introduced by Minkowski and Bouligand [48] to establish the dimension of an object using the measure theory and involves surroundings of each object point by additional border points with increasing size for the each iteration step. During the first dilation, each pixel is surrounded by a border of one pixel in width. Then, the reference element was a square of 3² pixels in size. At the second iteration step, each pixel was surrounded by a border of two pixels in width. The structuring element is then a square of 5² pixels in size [44]. At the third iteration step, each pixel is surrounded by a border of three pixels in width and the structuring element is then a square of 7² pixels in size.

In the third step, we measured the box-counting dimension of forested patterns after three dilations. With three dilations, all details referring to 1.078 km² areas were unclear, meaning that the smallest spaces disappeared. The pattern is fully connected in a fractal manner and the fractal dimension of the dilated forested area (D_Surf-dil3) should be higher than forest-free areas; D_Surf-dil3 should also be higher than 1. If D_Surf-dil3 differs substantially from D_Surf-dil0, the way of combining the smallest free spaces with the surrounding forest pixels follows a fractal organization very different from the spatial organization of the forest pattern at larger scales. On the contrary, if D_Surf-dil3 is not very different from D_Surf-dil0, the smallest free spaces do not follow a multi-scale pattern; they only influence the value of D_Surf-dil0 slightly by introducing noise into the measure [49].

The fourth step of the morphological characterization of the forest pattern was to measure the box-counting dimension of its dilated borders (D_Bord-dil3). The morphology of the borders is different from that of the areas and involved different ways to describe the geometry of forest areas. We used fractal analysis of dilated borders in order to eliminate the tiniest free spaces whose spatial forest organization is not well known and tends to be rather irregular.

Because of the forests’ fragmentation and heterogeneity, many borders of forest areas tend to exert irregularity with the D_Bord-dil3 greater than 1. For more compact patterns, D_Surf-dil3 tends to be higher than D_Bord-dil3 and vice versa, because the shapes of forested patterns are more diverse in the largest counties. This may lead to slight differences in the absolute values of the fractal dimensions, without affecting the general conclusions of the present paper.

3. Results

3.1. Spatial Evolution of Forest Covered Areas

Initially, the spatial evolution of forested, deforested and reforested areas was determined using ImageJ based on the Landsat GFC data (

Table 2. Characteristics of forest resources per development region during the period 2000–2012.

DRs/Parameter	W	Center	N-W	S-W Oltenia	N-E	S-E Dobrogea	S-Muntenia	Bucharest-Ilfov	Romania
FullDRs (km2)	32,028	34,100	34,159	29,212	36,850	35,762	34,489	1811	238,391
F2000 (km2)	13,468	14,719	12,327	9179	12,432	7061	7729	268	77,183
%F2000	42.1	43.2	36.1	31.4	33.7	19.7	22.4	14.8	32.4
F2012 (km2)	13,413	14,724	12,094	9182	12,122	7019	7617	266	76,438
%F2012	41.9	43.2	35.4	31.4	32.9	19.6	22.1	14.7	32.1
Def 2000–2012 (km2)	195	581	505	140	599	98	194	3	2315
Ref 2000–2012 (km2)	140	586	272	143	289	56	82	1	1570

FullDRs, Full areas of each DRs (km²); F2000, Forested areas 2000 (km²); %F2000, percent of forested areas; F2012, Forested areas 2012 (km²); %F2012, percent of forested areas; Def 2000–2012, Deforested areas 2000–2012 (km²); Ref 2000–2012, Reforested areas 2000–2012 (km²); W, West Development Region; Center, Center Development Region; N-W, North-West Development Region; S-W Oltenia, South-West Oltenia Development Region; N-E, North-East Development Region; S-E Dogrogea, South-East Dobrogea Development Region; S-Muntenia, South Muntenia Development Region; Bucharest-Ilfov, Bucharest-Ilfov Development Region;Data source: Data were extracted from images provided by [41].

). The forest-covered areas had generally decreased during the period 2000–2012. Only for a few DRs the observed reforestation was higher than deforestation (Center DR and S-W Oltenia DR).

3.2. Fractal Dimension of the Forested, Deforested and Reforested Areas (D_Surf-dil0)

At a national level, D_Surf-dil0 of forested areas was 1.50 which indicated an increased complexity, fragmentation and irregular arrangement by predominantly connected forested clusters. Due to the accentuated non-uniformity of the spatial distribution of the reforested and deforested areas, D_Surf-dil0 was as low as 0.93 for deforested and 0.88 for reforested areas respectively, indicating primarily detached forest clusters (

Table 3. Fractal dimension of forested, deforested and reforested areas (D_Surf-dil0) from development regions of Romania (2000–2012).

DR/FD	W	Center	N-W	S-W Oltenia	N-E	S-E Dobrogea	S-Muntenia	Bucharest-Ilfov	Romania
Forested areas
DSurf-dil0	1.63	1.62	1.58	1.58	1.54	1.44	1.43	1.14	1.50
Deforested areas
DSurf-dil0	1.05	1.15	1.08	0.98	0.99	0.83	0.86	0.53	0.93
Reforested areas
DSurf-dil0	1.00	1.16	1.05	0.98	0.91	0.77	0.80	0.40	0.88

Data source: Data were extracted from images provided by [41].

).

We observed that the West and Central development regions (DRs) had the highest values of D_Surf-dil0, due to the greater homogeneity of their forested areas, arranged in contiguous clusters, followed by North-West, West, South-West and North-East DRs. Scarcely forested DRs such as the South-East Dobrogea, South Muntenia and the Bucharest-Ilfov exerted a reduced values of D_Surf-dil0 which was lower than the national average. These DRs are counties characterized by limited areas of forest with a dominant relief of plains or low hills (Figure 3a). A clear difference was observed between the values of fractal dimension for forested areas as compared to deforested/reforested (Figure 3b).

The forests appear as small clusters, isolated patches or are distributed along the waters. Particularly, values were smaller (1.25) for the counties characterized by plains, as Brăila (BR) and Călărași (CL) and also for Bucharest (B) county (0.92) (Figure 4a).

The fractal dimension of deforestation indicates the fact that the highest values of D_Surf-dil0 are found in the Central, North-West and West regions which reveals a greater volume of exploited timber and that deforestation occurred in a rather unstructured manner. The lowest values of fractal dimensions from South Muntenia, South-East Dobrogea and Bucharest-Ilfov regions were caused by the limited areal extent of forests (Figure 4b).

The reforested areas display a smaller fractal dimension because these areas are smaller, but on the other hand they cover a large part of the fragmentation and irregularity induced by the deforestation. The fractal dimension D_Surf-dil0 of a reforested area is higher than the deforested area in South-West Oltenia and Center DRs where important reforestation campaigns were implemented. On the contrary, differences are observed between D_Surf-dil0 of deforested areas and reforested areas in the Bucharest-Ilfov region.

3.3. Dilated Fractal Dimension (D_dil3)

For the forested areas, D_dil3 provides interesting information about the fractal dimension obtained after dilation (Figure 5a,b). Two situations can be observed: the usual case with counties characterized by lower landscape (plains, low hills, sparsely forested counties) where the fractal dimension D_dil3 is less than D_Surf-dil0 (e.g., D_Surf-dil0 − D_dil3 for Galați (GL) was 0.32) and an atypical situation for counties with rugged landscapes (mountainous, hilly, densely forested) where D_dil was less than D_Surf-dil0 (e.g., D_dil3 − D_Surf-dil0 for Vrancea (VN) was 0.25) (Figure 6a). The main cause was that complex compact forestry areas markedly influence the D_dil computation.

This fact was obvious at the level of the following DRs: Bucharest-Ilfov, South-Muntenia with lower landscape and less forested areas revealing D_dil < D_Surf-dil0. The West and Center DRs with rugged landscapes, dominated by mountains or hills revealed D_dil > D_Surf-dil0. For deforested areas, the relation D_dil > D_Surf-dil0 disappeared and was replaced by D_dil < D_Surf-dil0. D_Surf-dil0 − D_dil3 varies between 0.2 in Bucharest (B) and 0.6 in Vâlcea (VL). At the DR level, the lowest values were found in Bucharest-Ilfov and South-Muntenia regions with the lowest extent of deforestations and in West and South-Oltenia DRs with the most deforestation. (Figure 6b). The same situation was observed in the reforested areas, where D_Surf-dil0 − D_dil3 varied between 0.18 in Bucharest (B) and 0.76 in Bihor (BH). The difference was larger than in the case of deforested areas due to less reforestation of forest areas (

Table 4. Fractal dimension of forested, deforested and reforested areas (D_dil3) from development regions of Romania between 2000 and 2012.

DR/FD	W	Center	N-W	S-W Oltenia	N-E	S-E Dobrogea	S-Muntenia	Bucharest-Ilfov	Romania
Forested areas
Ddil3	1.73	1.72	1.64	1.58	1.53	1.45	1.40	1.04	1.51
Deforested areas
Ddil3	0.56	0.70	0.67	0.46	0.56	0.41	0.46	0.17	0.50
Reforested areas
Ddil3	0.36	0.64	0.42	0.39	0.40	0.30	0.31	0.06	0.36

Data source: Data were extracted from images provided by [41].

).

At the level of DRs, the lowest values were recorded in South Muntenia, South East Dobrogea and Bucharest-Ilfov regions, where the reforestation was less pronounced, (while the highest values were recorded in West, North-West and South-West Oltenia DRs with significant reforestation. (Figure 6c).

3.4. Fractal Dimension of the Dilated Forest Areas (D_surf-dil3)

After applying a three-fold dilation (D_dil3) (Figure 7a), the fractal dimension D_Surf-dil3 of forest areas was increased compared to D_Surf-dil0 in all situations. The South-West Oltenia, North-East, North-West, Center, and West DRs with rugged landscapes, including mountains and high hills had a lower increase of D_Surf-dil3 compared to D_Surf-dil0 with a percentage increase with values between 9.6% and 23.9% according to the fractal analysis (Figure 7b).

For Bucharest-Ilfov, South-Muntenia, and South-East Dobrogea DRs with landscapes characterized by plains and low hills, there was a decrease of the fractal dimension by 4.5%–7.7%. The situation was much clearer at the county level (Figure 8a). The counties with mountainous landscapes and widespread forested areas (Harghita (HR), Suceava (SV), and Caraș-Severin (CS) counties) had lower D_Surf-dil3 values compared to D_Surf-dil0 (between 2.8% and 4.3%). In counties dominated by landscapes of plains/low hills with scarce forested areas (Bucharest (B), Călărași (CL), Teleorman (TR)), more significant increases of D_Surf-dil3 compared to D_Surf-dil0 were observed (with values between 11.3% and 36.6%).

It follows that as the difference between D_Surf-dil3 and D_Surf-dil0 increased, the fractal fragmentation also increased, especially for counties characterized by plain landscapes, because the artificial “fractalization” of space as the dilation was more pronounced.

At the level of DRs, D_Surf-dil3 of deforested areas increased proportionally with the degree of dilation, except for a higher increase of D_Surf-dil3 in South-West Oltenia towards the North-East DRs (

Table 5. The fractal dimension of forested, deforested and reforested areas (D_Surf-dil3) from development regions of Romania (2000–2012).

DR/FD	W	Center	N-W	S-W Oltenia	N-E	S-E Dobrogea	S-Muntenia	Bucharest-Ilfov	Romania
Forested areas
DSurf-dil3	1.70	1.69	1.68	1.69	1.66	1.57	1.58	1.39	1.62
Deforested areas
DSurf-dil3	1.42	1.50	1.42	1.36	1.35	1.22	1.25	0.99	1.32
Reforested areas
DSurf-dil3	1.43	1.50	1.42	1.36	1.34	1.22	1.25	0.99	1.31

Data source: Data were extracted from images provided by [41].

). This may be explained by the smaller areas of forest in Vaslui (VS), Iași (IS) and Botoșani (BT), mainly arranged in strips. At the county level, there was a very strong positive correlation between D_Surf-dil0 and D_Surf-dil3, with R² = 0.99, indicating that the deforested area was relatively small, compact and unevenly distributed.

Center and North-West DRs with varied landscapes, including mountains, hills and highlands had a lower increase of D_Surf-dil3 compared to D_Surf-dil0 (an increase of 31.1% to 32.6%). Bucharest-Ilfov, South-East Dobrogea, and South-Muntenia DRs with lower plains and hills showed a dilation induced increase of the fractal dimension with values between 47.5% and 109.8%. The mountain landscape of highly deforested counties (Suceava (SV), Vrancea (VN), Harghita (HR)) showed a lower increase of D_Surf-dil3 compared to D_Surf-dil0 with values between 19.6% and 34% (Figure 8b). In counties dominated by the lower plains and hills landscapes, and with a low degree of forest areas such as Bucharest (B), Constanța (CT), and Vaslui (VS) the increases of D_Surf-dil3 compared to D_Surf-dil0, were larger with values between 54.3% and 161.5%. D_Surf-dil3 of reforested areas increased proportionally to the dilation, with the exception of the higher increase in the West DR compared to the North-West DR (due to higher values of D_Surf-dil3 in Caraș-Severin (CS) and Timiș (TM)), because the reforestation has been conducted in small and relatively uniform patches. At the county level, there was a strong correlation between D_Surf-dil0 and D_Surf-dil3 (R² = 0.83), but this correlation was weaker for deforested areas because the reforestation was distributed more chaotically. Regions with varied landscapes, including mountains, hills and highlands such as South-West Oltenia, North-West and Center showed a lower increase of D_Surf-dil3 compared to D_Surf-dil0 between 30.5% and 40.2% compared to DRs with landscapes of lower plains and hills such as Bucharest-Ilfov, South-East Dobrogea, and South Muntenia, where the dilation induced increase of fractal dimensions ranged between 58.7% and 211.8%. The values of D_Surf-dil3 − D_Surf-dil0 for reforested areas were much higher than the D_Surf-dil3 − D_Surf-dil0 for the deforested areas because of the high extent of fragmentation. However, the counties with landscapes dominated by mountains and dense forests (Gorj (GJ), Suceava (SV), Maramureș (MM)) showed a lower increase of D_Surf-dil3 compared to D_Surf-dil0 (from 35% to 18.3%). In counties with landscapes dominated by lower plains and hills with limited areas of forest (Bucharest (B), Ilfov (IF)), a more pronounced increase of D_Surf-dil3 was detected compared to D_Surf-dil0 with values between 61.3% and 334.1% (Figure 8c).

3.5. Fractal Dimension of the Dilated Forest Areas Borders (D_Bord-dil3)

Extracting the border of three-fold dilated images _dil3 allowed for the calculation of the fractal dimension D_Bord-dil3. D_Bord-dil3 provided information about the fractal shape of forested area borders and it increases when was growing the irregularity of the borders (Figure 9 a,b).

At the county level, the smallest fractal dimensions of forested areas (indicating lower fragmentation and complexity of the border) were found for forested areas of Bucharest (B), Vrancea (VN) and Caraș-Severin (CS) (very small but compact) whereas the largest fractal dimensions were found in Vaslui (VS), Galați (GL), Iași (IS) as the forested areas were distributed as elongated strips (Figure 10a).

Lowest values of D_Bord-dil3 were recorded in South-East Dobrogea and Bucharest-Ilfov where deforestation was sporadic and relatively compact, the highest in N-E where the fragmentation of forest areas was intense due to deforestations and South-West Oltenia. At the county level, for Bucharest, but also for Vrancea (VN), Covasna (CV), and Caraş-Severin (CS) counties for deforested areas, the lowest D_Bord-dil3 indicated fragmentation and less complexity due to very small but compact forested areas (where the deforestation was of a compact shape) and the highest fractal dimensions were found in Vaslui (VS), Cluj (CJ), and Iași (IS) counties (where the deforestation was shaped as elongated strips) (Figure 10b).

At the level of DRs, the lowest values of D_Bord-dil3 were recorded in South-Muntenia, South-East Dobrogea and Bucharest-Ilfov regions, with small and relatively homogeneous deforested areas. The highest values of D_Bord-dil3 were registered in DRs with widespread deforestation of more irregular shape due to the configuration topography, and with forest poaching, such as the Center, West and North-West. The lowest values of D_Bord-dil3 of reforested areas were calculated, similarly to deforested areas, in South Muntenia, South-East Dobrogea and Bucharest-Ilfov DRs, with low and relatively homogeneous reforestation, especially by plantation. The highest values of D_Bord-dil3 were registered in DRs with fragmented landscapes (e.g., Center, West and North-West DRs), where self-regeneration and plantation campaigns were imposed (

Table 6. Fractal dimension of forested, deforested and reforested areas (D_Bord-dil3) from development regions of Romania (2000–2012).

DR/FD	W	Center	N-W	S-W Oltenia	N-E	S-E Dobrogea	S-Muntenia	Bucharest-Ilfov	Romania
Forested areas
DBord-dil3	1.29	1.29	1.32	1.33	1.35	1.28	1.31	1.22	1.30
Deforested areas
DBord-dil3	1.30	1.35	1.30	1.26	1.24	1.13	1.15	0.91	1.21
Reforested areas
DBord-dil3	1.32	1.33	1.31	1.28	1.21	1.13	1.18	0.84	1.20

Data source: Data were extracted from images provided by [41].

).

At the county level, the lowest fractal dimension values (between 0.57 and 1.14) were calculated the borders of forest areas in the regions of plains or low hills. e.g., Bucharest (B), Teleorman (TR), Botoşani (BT), and Vaslui (VS) counties, indicating their fragmentation and lower complexity. In contrast, the highest fractal dimension values of 1.36 and 1.44 were calculated for Alba (AB), Argeș (AG), Harghita (HR), and Suceava (SV) counties (Figure 9c).

4. Discussion

4.1. Forest Changes in Romania

Our results showed important spatial forest changes in Romanian DRs during recent decades. By using a satellite based information on forest changes and fractal analysis, we provide information describing the state and changes in forest areas in all eight DRs of Romania during the period 2000–2012. The satellite based per-pixel assessment of forested, deforested and reforested areas was complemented by fractal analysis revealing morpho-structural and textural differentiations of forest disturbance in different DRs, providing a quantification of changes in the uniformity, homogeneity/heterogeneity and compaction/fragmentation of forests.

Forest areas in Romania are characterized by a general downward trend over the period covered by the study due to an increased level of legal and illegal logging [50]. This development is reflected by the economic changes during this period and is spawned by changes in legislation which encouraged deforestations [11,36,37,38]. In a longer historical perspective, forest cover in Romania has been characterized by significant changes over the last 100 years, mostly due to socio-political factors that have had significant impact on the forest management [51]. The 1926 to 1948 period has been characterized by a decrease in private land ownership, whereas from 1948 to 1989, in the socialist period, the forests were entirely state-owned. In the post-socialist period, an increase in private lands is a result of the three privatization laws from 1990 to 2004 [27,52]. Logging was particularly intense during the Socialist-era in the 1960s and 1970s but forest harvesting declined during the 1990s as a consequence of the collapse of the timber market and lower institutional strength and weak law enforcement, triggering an increase in illegal logging at the expense of planned logging [27,51,52]. Areas of primary forest in Romania covered aproximately 4000 km² in 1984 but this was reduced to 2185 km² in 2004 [28,53] and consequently the habitats for species were also diminished [24].

4.2. Usage of Fractal Analysis for Quantification of Forest Changes

We have evaluated the forest resources [41] using the box-counting method [44], determining the fractal dimension of the forest areas (D_Surf-dil0), the fractal dimension of the dilated forest areas (D_Surf-dil3), the fractal dimension of the border of the dilated forest areas (D_Bord-dil3) and dilation fractal dimension (D_dil3). Our results add new insights to forest cover change as compared to previous studies based on fractal analysis [3,11,13,36,37,38]. The findings of this study clearly confirmed the hypothesis that fractal analysis can provide an interesting new way of quantifying the degree of fragmentation and organization of forest areas characterized by deforestation and reforestation. Here, we used fractal analysis to study the characteristics of the forestry morphology for various levels of forested, deforested and reforested areas in Romania. Fractal analysis was used because of the ability of the fractal dimension to provide a quantitative description of the forest environment. By measuring the contrast of forestry patterns across scales, this method is applicable for the evaluation and characterization of forest area types. Through fractal analysis, we obtained complementary information on deforestation and reforestation/regeneration compared to classical change detection analysis based on image classification because the fractal analysis is able to analyze irregular spatial structures.

In the fractal analysis of forested areas, we used the box-counting method to describe spatial patterns [11,36,37,38] and the dynamics of deforestation [34]. In previous studies, fractal analysis of forest areas was determined by using the Fractal Fragmentation Index [13,36] or Minkowski Dimension and Pyramid Dimension [33]. Until now, the dilation method was used to evaluate the quality of the commune or urban built-up areas [44] and we clarified the importance of utilizing the dilation method in understanding the textural organization of forest areas and the effects generated by deforestation and reforestation. The main findings are summarized and discussed in the following five points:

• We show that fractal analysis can offer valuable textural information about the spatial organizations of forest resources. Also, the fractal analysis highlights the spatial effects of deforestation on forest resources.
• D_dil3 was shown to be a useful index that indicates the degree of compactness of forested areas. D_dil3 showed higher values than D_surf-dil0 only in the situation of densely forested areas and with a great degree of compactness. Otherwise, for areas less forested, deforested or reforested, D_dil3 was less than D_surf-dil0.
• The difference between D_surf-dil3 and D_surf-dil0 can be utilized as an index of the forest area complexity organization. This indicates that the value of D_surf-dil3 − D_surf-dil0 increases for more heterogeneous, fragmented and scattered forest areas. If forest areas were compact and organized in big clusters, the difference between D_surf-dil3 and D_surf-dil0 tended to be lower.
• D_Bord-dil3 offered important textural information about deforestation and reforestation. The more compact the deforestation or the reforestation, the lower was the value of D_Bord-dil3.
• The fractal analysis highlighted morpho-structural and textural differences between forested, deforested and reforested areas grouped by DRs or counties. Differences were observed between landscapes dominated by mountains and high hills (more forested and compact organization) as compared to landscapes of lower plains or hills (less forested, more fragmented, with many small and isolated clusters). Also, obvious differences were found between morpho-structural and textural patterns of forested areas (more compact) as a function of deforested and reforested areas (more isolated points, patches randomly distributed). At the level of DRs, areas of deforestation not followed by a sufficient reforestation show a change from a compact structure of forest areas to a chaotic structure with multiple free spaces that are non-forested pixels.

4.3. Perspectives for Sustainable Forest Management

The current study shows the importance of utilizing the dilation method in understanding the textural organization of forest areas and the effects generated by deforestation and reforestation. The method, based on fractal analysis, thereby shows the advantage of revealing the effects and patterns of forest disturbance altogether. Fractal analysis of both D_Surf and D_bord is involved in characterizing forestry morphology.

Modelling of deforestation through fractal analysis, targeting forest textural differentiation, provides interesting complementary information to traditional per-pixel deforestation mapping. By combining such methods, it is expected that decision makers can gain an improved starting point for understanding the underlying driving forces of deforestation, which is instrumental for optimal planning of sustainable forestry and future exploitation of forest resources. Specifically, the morpho-structural and textural information on forest fragmentation and heterogeneity provided by fractal analysis can be of added value in assessing the risk of landslides, floods and flash floods, which is higher in the watersheds from DRs characterized by high deforestation rates in the Carpathians and Subcarphatians [21,22]. Also, forest areas have an important ecological function in relation to ecosystem services such as biodiversity. Here, increased landscape fragmentation is known to be a major determinant for a biodiversity decline [54,55] and fractal analysis offers a powerful way of quantifying landscape or forest fragmentation and assessment of temporal changes herein from remote sensing time-series.

However, the model did not sufficiently highlight the complexity of morphology of forest areas at a local scale, such as the ratio of clusters and isolated clumps of forest. This shortcoming may be resolved by a complementary model of local fractal analysis including mass-dimension (for areas) and ruler-dimension (for the border of the principal cluster).

It must be kept in mind that the Landsat based global forest cover (GFC) maps [41] used as an input for the fractal analysis are based on definitions of forest loss. Forest loss in the GFC data is defined as the disturbance or complete removal of tree cover canopy (below 25% tree canopy cover) including any conversion of natural forests by local communities, also including deforestation for establishment of plantations. Therefore, the forest mapping should be intersected with ground knowledge for the correct interpretation of patterns of forest loss. Moreover, the Landsat based forest cover maps are based on optical remote sensing that is known to suffer from shortcomings in the mapping of forest degradation, mainly due to the lack of sensitivity of optical instruments to changes within and below the forest canopy [56].

Finally, a more detailed analysis, e.g., differentiating the type of forests (deciduous or coniferous forest, private or state forests) might help in establishing of the novel forest management which would be better suited for practical implementation.

5. Conclusions

Fractal analysis of forest areas provides the possibility to calculate the fractal dimension both of shapes and areas. Therefore, the fractal analysis is a tool to evaluate spatially explicit geographic phenomena in more detail and can improve our knowledge of the spatial organization of disturbances in forested areas. We analyzed changes in forest resources from Romanian development regions affected by both deforestation and reforestation by using fractal analysis. Through calculation of four forest area fractal dimensions we generated the quantitative morpho-structural and textural information about the degree of uniformity, fragmentation, heterogeneity and homogeneity complementary to traditional non-textural estimates of changes in the spatial extent of deforestation and reforestation. Therefore, we recommend the application of fractal analysis to the time-series of satellite remote sensing data as a new approach to providing additional information for the improvement of forest area management.