*Article* **Determining Long-Term Land Cover Dynamics in the South Baltic Coastal Zone from Historical Aerial Photographs**

**Andrzej Giza 1,\* , Paweł Terefenko <sup>1</sup> , Tomasz Komorowski <sup>2</sup> and Paweł Czapli ´nski <sup>3</sup>**


71-101 Szczecin, Poland; pawel.czaplinski@usz.edu.pl

**\*** Correspondence: andrzej.giza@usz.edu.pl; Tel.: +48-91-444-24-32

**Abstract:** Coastal regions are dynamic environments that have been the main settlement destinations for human society development for centuries. Development by humans and environmental changes have resulted in intensive land cover transformation. However, detailed spatiotemporal analyses of such changes in the Polish Baltic coastal zone have not been given sufficient attention. The aim of the presented work is to fill this gap and, moreover, present a method for assessing indicators of changes in a coastal dune environment that could be an alternative for widely used morphological line indicators. To fulfill the main aim, spatial and temporal variations in the dune areas of the Pomeranian Bay coast (South Baltic Sea) were quantified using remote sensing data from the years 1938–2017, supervised classification, and a geographic information system post-classification change detection technique. Finally, a novel quantitative approach for coastal areas containing both sea and land surface sections was developed. The analysis revealed that for accumulative areas, a decrease in the land area occupied by water was typical, along with an increase in the surface area not covered by vegetation and a growth in the surface area occupied by vegetation. Furthermore, stabilized shores were subject to significant changes in tree cover area mainly at the expense of grass-covered terrains and simultaneous slight changes in the surface area occupied by water and the areas free of vegetation. The statistical analysis revealed six groups of characteristic shore evolutionary trends, of which three exhibited an erosive nature of changes. The methodology developed herein helps discover new possibilities for defining coastal zone dynamics and can be used as an alternative solution to methods only resorting to cross sections and line indicators. These results constitute an important step toward developing a predictive model of coastal land cover changes.

**Keywords:** land cover; dune coast; air photograph; South Baltic Sea

#### **1. Introduction**

Coastal zones are highly dynamic regions that exhibit unique atmospheric, hydrospheric, lithospheric, and biospheric characteristics. Numerous studies have indicated that weather conditions and rising sea levels are the primary factors influencing coastal erosion [1].

Since the mid-twentieth century, increased cyclonal activity has been observed in winter periods in the North Atlantic due to global climate change [2,3]. The number of extreme storm surges in the Baltic Sea is increasing steadily [4,5], which along with milder winters and limited ice cover, further exacerbates coastal erosion [6]. The observed rates of change in coastlines, particularly in case of erosion, have been a significant concern for communities inhabiting littoral zones. Appropriate management closely linked to a balanced development of littoral areas is becoming a major challenge for the research community owing to its role in the knowledge-based shaping of environmental policies [7]. The study of long-term changes in land cover is important, since it supplements the knowledge of

**Citation:** Giza, A.; Terefenko, P.; Komorowski, T.; Czapli ´nski, P. Determining Long-Term Land Cover Dynamics in the South Baltic Coastal Zone from Historical Aerial Photographs. *Remote Sens.* **2021**, *13*, 1068. https://doi.org/10.3390/ rs13061068

Academic Editor: Feng Ling

Received: 10 February 2021 Accepted: 10 March 2021 Published: 11 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 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 (https:// creativecommons.org/licenses/by/ 4.0/).

coastal zone dynamics by analyzing vegetation cover, which is a crucial element stabilizing shore ecosystems. Changes occurring in the coastal zone landscape and its ecosystems are critical in terms of spatial planning and determining natural environment conditions.

In recent years, remote sensing (RS) and Geographic Information Systems (GIS) have gained in popularity, in assessing the spatial and temporal dynamics of coastal zones, as powerful and cost-effective tools [8–10]. RS methods have been particularly useful in analyzing land cover changes. Satellite images constitute a basic source of data, providing information on environmental conditions over a far greater area than that acquired by aerial imaging. In turn, aerial photographs are a valuable source of information regarding environmental conditions at a given moment with a high degree of accuracy. RS methods and spatial information systems are valuable when documenting and measuring changes to the terrain. Adopting a dynamic perspective toward the landscape promotes mutual human and nature interactions [11]. Furthermore, remote sensing tools are successfully applied in studies for detection changes in shallow nearshore areas [12–14]. The use of archival photos taken using monochromatic technology offers many possibilities, as confirmed by numerous studies [15–21]. This technology requires from the photograph interpreter significant experience and knowledge of the physical and geographical conditions of the analyzed area. Aerial photographs are also employed for monitoring changes in seashore dynamics; however, such studies are conducted by analyzing the location of the dune or cliff base line in individual years. Currently, precise measurements are commonly made using satellite navigation, taking advantage of real time kinematic GPS to generate transverse profiles of the coastal zones. For larger areas, light detection and ranging (known as LiDAR) technology has been successfully employed [22–25].

The South Baltic land cover changes and their dynamics have long been a matter of widespread research interest. The first detailed study of the South Baltic coast concerned the Island of Wolin's flora [26]. The paper listed species and their locations, as well as plant communities of many different species, both common and rare. Short and Hesp [27] created a classification of dune shores based on the percentage of vegetation cover, in which the morphological development was determined by the vegetation affecting eolian processes and form creation [28]. Coast development on the tideless Baltic Sea shores, and the possibilities for forecasting in the area, were also explored by Musielak and Furma ´nczyk [29]. Furthermore, Dudzi ´nska-Nowak analyzed dune/cliff base changes between 1938–1996 within the Rewal Commune [30], while the studies conducted by Łabuz regarding the impact of anthropogenic activity on dune vegetation cover on the Island of Wolin between Mi ˛edzyzdroje and Swinouj´scie confirmed the natural trends ´ of vegetation cover development in that area [31,32]. Thus far, most studies have been time- or site-specific. Spatiotemporal changes in land cover have not been sufficiently addressed. The most extensive study to date has been presented by Bielecka et al. [33]. Her paper provides the first comprehensive quantitative and qualitative analysis of land cover changes and their characteristics in the Polish Baltic coastal zone over the last three decades. It does not, however, offer local-scale information, since the analysis' scale is adjusted to the resolution of the European CORINE Land Cover database.

Therefore, to fill this gap the present study offers analyses of the coastal dynamics on a local scale, focusing on long-term land cover changes. A new reproducible quantitative approach using basic field analysis to define the evolutionary trends of individual sections of a dune coast is proposed. For this purpose, a spatial database of archival aerial photos was developed and subsequently used to analyze changes in land cover. Finally, models of seashore evolutionary trends were created based on the results of the land cover condition analysis. These methodologies constitute an important step toward developing a predictive model of coastal land cover changes. Furthermore, the implementation of this new approach will enable researchers to test scientific hypotheses stating that the analysis of land cover changes in enumeration units could be an alternative for widely used morphological line indicators in a coastal dune environment.

#### **2. Materials**

#### *2.1. Study Site*

The Baltic Sea is a non-tidal, semi-enclosed, and shallow body of brackish water with a coastline divided between 14 countries. Poland has a 500-kilometer-long diversified coast [34]. In several locations along the coast there are sections of cliffs that account for 18% of the entire coastline [5]. Apart from small alluvial sections (less than 3% of total), the remainder is a spit- and barrier-type coast, with dunes ranging in height from less than 2 to 49 m.

Generally, those dunes are subject to erosion along the whole Polish coast, and thus usually long sections of the dune coast have to be reinforced with some technical protective infrastructure. Several locations that are exceptions to this, such as the Spit of the Swina ´ River or Swinouj´scie region (both in the western part of the Polish coast), experience ´ shore accumulation. These unique sites were chosen as the subject-matter of our land cover change analysis. The research area encompasses the western part of the Polish coast with several regions of dune coast located within the administrative limits of the West Pomeranian Voivodeship. It comprises a stretch from the city of Swinouj´scie to Mrze ´ zyno ˙ sea resort (Figure 1). To ensure accuracy and verifiability, and to allow for a comparison of our results with those of previous studies, we narrowed the scope of this study to the dune-covered coast. Implementing such limitation rules resulted in several exclusion areas of varying sizes.

**Figure 1.** Location of the investigated area.

The first excluded fragments represent coastal cliffs. The decision to drop these parts of the coast was due to the substantial differences in height occurring between the cliff edge and base. Accounting for these differences would require the application of different analytic techniques to determine land cover changes.

Furthermore, considering that the analyzed terrain comprised areas where natural processes of shore and land cover evolution occurred, the analysis also excluded all coastal segments comprising river estuaries (the mouths of the Swina and the Dziwna rivers) ´ with their breakwaters, as well as areas modified by intensive construction. Moreover, the exclusion group comprised areas with strong shore fortifications, heavy-duty concrete

bands, shore reinforcements with artificial dunes and plantings, and places featuring anthropogenic changes located in the direct vicinity of the dune. Finally, human-modified spaces such as military grounds were also excluded from the analysis. As substantial changes in the area concerned were caused by military activity, they were masked in the 1973 aerial photograph that was available to us, thus preventing comparisons with other photos that we had.

The climatological conditions characterizing the investigated area are controlled by atmospheric pressure patterns and highly influenced by the Atlantic Ocean. The wave energy in the investigated area is medium, with an average significant wave height calculated for the 1958–2001 period ranging from 0.37 m in summer to 1.04 m during winter months. Offshore waves are mainly driven by westerly and southwesterly winds [35].

In recent decades, the highest absolute amplitude of sea level changes in the study area was recorded in 1984 and reached a value of 2.79 m, while the most extreme storm surge (a combination of the sea level and the waves) occurred in November 1995, with a water level of +1.61 m above mean. However, extreme value analyses showed that a 100-year storm surge in the western part of the Polish coast could reach +1.71 m above mean, and a 500-year event would exceed 2 m [5].

The coastal zone itself appears not to be crucial to the national economy, as the maritime sectors provide employment to less than 1% of the total working population. Additionally, the region is rather scarcely populated. The population of the municipalities enjoying direct access to the sea accounts for less than 3% of Poland's entire population. The average population density in the municipalities in the investigated areas is a maximum of 95 persons per km<sup>2</sup> , far below the national average of 123 [33]. However, the coast is by far the most popular vacation destination in Poland. It is estimated that the coastal zone attracts approx. 35% of all tourist traffic in Poland during July-August. This corresponds to an average of 71 tourists per 100 residents during the summer [5]. Furthermore, the Baltic Sea coastal zone boasts more than a 20% share in the country's high nature value areas, which attract tourism development and improve the residents' quality of life [33].

#### *2.2. Data*

In this study, five series of RS images from 1938, 1951, 1973, 1996, and 2017 were used. The 1938 photomap is in the worst condition due to its age. Although it contains a map graticule along with location names and object descriptions, those elements had a limited impact on the image interpretation outcome. Archival panchromatic aerial photographs from 1951 and 1973 also exhibit signs of damage, such as scratches and discoloration, which directly affects image clarity, although the photographs from 1951 offer satisfactory brightness and clarity, while those from 1973, despite their smaller scale, provide a better contrast. The photographs taken in 1996 are in the best technical condition, as they bear no signs of use and are the most suitable for interpretation. Finally, digital color orthophotomaps from 2017 were obtained from the Head Office of Geodesy and Cartography.

All the 117 aerial photographs analyzed had similar scales and pixel sizes, constituting a valuable comparative material for the detailed analysis of land cover changes. The numbers of images per series, as well as their scales, representing the spatial range of the image, their scan resolutions, and the resulting spatial resolutions (pixel sizes) of the digital versions are presented in Table 1.


**Table 1.** Archival aerial photograph details.

#### **3. Methodology**

The methodology comprised three basic steps: image processing, land cover detection, and statistical analysis of the changes. The methodology workflow is presented in Figure 2.

**Figure 2.** Methodology workflow.

#### *3.1. Image Processing*

A photomap from 1938 and archival aerial photographs from 1951 and 1973 were subjected to geometric correction (rectification) into a system of coordinates of the orthophotomap compiled with the use of color photographs taken within the scope of the PHARE program in 1996, which organized all the data in a uniform Pa ´nstwowy Układ Współrz ˛ednych Geodezyjnych 1992—Polish National Coordinate System 1992 (PUWG 1992) layout. The process involved combining each image with a system of coordinates using n-th degree polynomials [36,37]. Rectification was performed using first degree polynomial transformation in Erdas Imagine 8.5. The purpose of the method is to ensure that an unprocessed image can provide a suitable representation.

The parameters defining a mathematical equation between the previous and the new image reference frame are defined using the ground control points (GCPs), which are points on an output image that correlate to points on the input image to facilitate transformation. The minimum number of ground control points required for the first-degree polynomial is three [38]. However, during rectification, many GCPs were marked and distributed evenly over the image surface, including the marginal areas. Overall, 143 points were placed over the entire surface. The process was performed for each photograph series, beginning with the most recent photos, i.e., those taken in 1973, followed by those dating from 1951 and 1938.

Image rectification is crucial due to the impact of the cartometricity of the processed image on the subsequent analyses. As a result of rectification, the accepted root-meansquare error (RMSE) value in the case of the available material ranged between 2 and 5 m, while the greatest errors were observed for the photomap of 1938. This was so because the image provided for our analysis was processed into a digital form with the use of a flat scanner with A3 formatting. Consequently, the original photomap sheets were scanned in two parts, dividing each sheet into the E (eastern) part and the W (western) part. The RMSE values for the remaining photographs from 1951 and 1973 were within the range of 2 to 3 m.

#### *3.2. Land Cover Detection*

For a detailed analysis of land cover changes, an approx. 50-km-long area of dune sections of the coastline was designated as shown in five photograph series spanning the following periods: 1938–1951, 1951–1973, 1973–1996, and 2017–1996. ArcGIS Pro 2.6 software was used, and the entire area of interest was divided into basic fields (Figures 1 and 3). Each such subarea represented a 1-km-long and 300-m-wide coastal fragment. On several occasions, basic field lengths were slightly shortened to border areas excluded from the analysis. In total, 51 basic fields were extracted and analyzed (Figure 3).

**Figure 3.** An example of land cover detection for basic field no. 9.

Aerial photograph interpretation was based on Stone [39], who discussed the images of individual plant species captured in photos. How vegetation cover is represented in aerial photographs depends on a multitude of factors, which hinders interpretation and limits the possibility of using standard keys. Vegetation cover constitutes an environmental component that is typically clearly visible, even in panchromatic aerial photographs. However, it is not always distinctly interpretable or identifiable [40].

Within the scope of the supervised classification of the archival aerial photograph series, various states of land cover divided into classes were obtained separately for each analyzed area. Firstly, basic fields were divided into two smaller fragments, isolating an area comprising the sea and the beach, as well as another area comprising the beach, dune vegetation, trees, and shrubs. This particular division stemmed from the fact that in 8-bit photographs, the pixel phototones for water-covered and vegetation-covered areas are remarkably similar. Therefore, to ensure the correct interpretation and classification of archival photograph details, the separate class of "water" had to be singled out.

Once the class of "water" was distinguished, the remaining objects were classified in terms of recognizable properties. Direct object properties were considered, such as image shape and size, phototone, color, texture, and structure. Indirect characteristics were also important, as they indicated the presence of objects or their properties, such as own shadow, shadow cast, topographic location, and pixel values related to other landscape elements, which were verified on the basis of photograph image structure and texture [41].

Identification features may also be described using tools provided by image processing systems. In a digital image, texture is understood as the superficial pixel distribution of various degrees of brightness values, accounting for the regularities and mutual relationships occurring in that distribution [42]. In object recognition, the influence of an interpreter's subjective assessment may be significantly decreased. A digital image record allows for the precise definition of the nature of tone changeability and for the capture of subtle differences in the reflection of the individual ranges of electromagnetic radiation [41].

A supervised classification was done to locate specific areas within the image that show homogenous categories of the land cover types. Object base classification was used with a manually created training sample for each image over identifiable cover types. The process involved using the Training Samples Manager tools in the ArcGIS PRO software in order to set the neighborhood mode, the geographic constrains of the pixels, and the spectral Euclidean distance. Multiple areas were drawn and grouped for a single category. Those areas represented the best match of pixels. The process was repeated for every known category on the image. The supervised classification used the maximum likelihood parametric classification rule.

When interpreting vegetation cover, we were able to differentiate between two basic groups: tree vegetation (with shrubs) and grasses. Non-forest formations can be discerned through changes in texture, structure, tone, and color to determine the different species; the density of vegetation cover itself is essential as well. Grass vegetation is typically smooth and nearly uniform in color, while coastal grasses are lighter than other vegetation. Eventually, the following land cover classes were specified:


These classes have inclusive names to account for a wide range of species. The term "grass" indicates low-growing vegetation up to a height of approximately 0.5 m, comprising various plant communities. The term "tree" indicates vegetation growing higher than 0.5 m and contains shrubs and young trees. Both these classes were clearly distinguishable in our aerial photographs.

This stage was also important for collecting the testing sites of the study area. These are very important to determine the accuracy assessment for each classification algorithm and to check the validation of classifications. Usually, whenever possible, the testing sites should be represented as ground control points collected from the field of the study area. Because the presented work is based on the analysis of historical data, such control points (one for each of 4 classes) has been collected using office work from each of the 117 analyzed photos.

The final step of the image classification was the accuracy assessment stage. This process was performed as an estimation with the aid of a remotely sensed dataset for classification conditions, and it is useful for the evaluation of the classification approach. It is also important to determine the error that might be involved. Within the presented work, a confusion matrix approach has been implemented [43]. Control points collected randomly for each class on each image have been set against the achieved results, thus enabling the validation of the results of our classification. The overall average accuracy of classification was 89.5%. The values were different for each series of photographs, reaching 82.43% for 1938, 90.25% for 1951, 92.22% for 1996, and 93.12% for 2007, respectively.

#### *3.3. Statistical Analysis*

Having identified and measured the area covered by the individual classes, the results were recalculated from real values expressed in square meters to relative values expressed as a percentage share of each class. These shares allowed for a detailed analysis of land cover changes without comparing the sizes of basic field surfaces in relation to their actual surfaces. Consequently, a change inferred from the share size of a given class could be analyzed both for an individual basic field and the entire analyzed section of the coast extending from Swinouj´scie to Mrze ´ zyno. ˙

All statistical calculations and analyses were realized within Statistica 12.5 software. Firstly, in order examine changes over time, these land cover classes were subjected to statistical analysis. First, several standard descriptive statistics, histograms, and standard box plots were used to indicate variability in the land cover classes. The box plot, also called box-whiskers plot, was used in order to present the variability of the land cover classes. Values larger than the sum of the third quartile and three quarter deviations (Q\_3+3 (Q\_3- Q\_1)/2), as well as those less than the difference of the first quartile and three quarter deviations (Q\_1+3 (Q\_3- Q\_1)/2), are considered to be atypical, i.e., outliers [44].

Further on, the share values of the individual land cover classes were subjected to cluster analysis, which is used to arrange objects in groups according to the degree of relation [44]. This method is most helpful in the exploratory phase of research, where no hypothesis has yet been formulated. Therefore, cluster analysis does not constitute a statistical test but helps group objects according to their similarities. In the present research, it was used to cluster two different types of data. First of all, it was employed to distinguish groups of basic fields with similar shares of land cover classes. Second, it was used to examine and create a clustering process of similar patterns in land cover changes between the time periods studied.

To this effect, the Euclidean distance—the geometric distance in a multi-dimensional space—was used to determine the distance between the new clusters. The Ward's method that was used also assumes that every object constitutes a separate cluster. At every subsequent stage, groups are created by combining objects and clusters formed at earlier stages. Consequently, the objects were combined into ever-growing clusters until a sufficient level of grouping was achieved. The method employs variance analysis to estimate the distance and aims to minimize the sum of the deviation squares inside clusters. At every stage, of every possible pair of clusters that can be combined, the cluster that eventually yields the minimum diversification is selected [45]. The measure of such diversification in relation to average values is the error sum of squares, which is defined by the following formula:

$$ESS = \sum\_{i=1}^{k} (x\_i - \ x)^2 \tag{1}$$

where:

*x*—the base value of the variable constituting the segmentation criterion,

*xi*—value of the variable constituting the segmentation criterion for the i-th object,

k—number of objects in a cluster.

Linear regression, being the simplest regression variant, was employed to determine the development trends of individual classes of land cover over time. It assumes that there exists a linear relationship between an explained variable and an explaining variable. In linear regression, it is presumed that an increase of one variable is accompanied by the growth or decrease of the other variable. A regression function takes on the form of a linear function:

$$\mathbf{y} = \mathbf{b}\mathbf{x} + \mathbf{a} \tag{2}$$

where:

a and b are regression estimators.

Linear regression analysis is used to determine such coefficients with which the model could best predict the value of a dependent variable. Such approach ensures that the estimation error is as small as possible. It also enables us to set a development trend of land cover class changes which could be used for defining model data in further analyses [44].

#### **4. Results**

#### *4.1. Descriptive Analysis of Land Cover*

The percentage shares of the individual classes are presented in the form of a stacked column chart (Figure 4). Characteristic arrangements were observed for areas 1 to 15, which comprised the section stretching from Poland's western border through the mouth of the Swina River to Mi ˛edzyzdroje. This section exhibited an accumulative character with ´ a fragment of a stabilized shore lying to the west of Mi ˛edzyzdroje (the areas 13 to 15). In the aerial photographs from 1938 and 1951, areas 1 to 12 revealed a wide beach along with grass vegetation covering the territory; there were no visible trees or shrubs there, though. Trees and shrubs in those areas appeared as late as 1973, and their surface area expanded greatly in 1996. This territory clearly implied accumulative processes, occurring particularly in areas 5, 6, and 7, lying to the east of the breakwater at the mouth of the Swina River. In area 5, water occupied 66.2% in 1938 and only 36.3% in 1996, whereas in ´ area 6 water area shrank from 53.9% in 1938 to 28.3% in 1996 (Table S1).

The areas between 11 and 15 constituted a section of a stabilized shore. The accumulation found in that stretch of the coast was slight, but a natural succession of vegetation cover entering new areas was observed. In areas 13 and 14, trees and shrubs were recorded in 1938 and 1973. The shore was stabilized, and the water surface area was virtually unchanged, accounting for approximately 60%. The proportions of the remaining land cover classes kept changing. An expansion of tree and grass cover at the expense of surface clear of vegetation could be observed. Areas 16 to 37 were clearly more stable in terms of land cover change dynamics (Tables S2 and S3). Between areas 16 and 23 (Tables S1 and S2), lying to the west of the mouth of the River Dziwna, the water surface area exhibited minimal to no expansion, while changes in the land cover structure were noted. Such changes were particularly visible between the periods 1938–1951 and 1951–2017. During the former, the surface clear of vegetation expanded along with an increase in the beach width. During the latter, a rise in grass vegetation became evident, followed by the appearance of trees and shrubs, which were visible in those areas in 1996.

The portion of the coast located between Dziwnów and Trz ˛esacz, comprising the areas from 24 to 37, was a stabilized shore fragment. Water and beach surface areas remained the same, and the tree area increase became noticeable with a simultaneously shrinking surface occupied by grass.

In areas 38 to 51 (Pogorzelica-Mrzezyno), the water surface area advanced, which ˙ implicated that erosion processes were underway (Table S3). Significant areas clear of vegetation were observed, specifically in 1938 and 1951 (Figure 4). In subsequent years, a rise in the surface occupied by trees, shrubs, and grasses was observed.

The investigated status of land cover, with account taken of particular classes within the analyzed coastal sections, revealed that the terrains extending in the direct vicinity shared similar proportions of land cover structure. The dynamics of the changes occurring in those terrains were also similar. Additionally, it was noted that the territories excluded from the research, such as river estuaries or cliff-coast sections, were naturally distinct areas with different land cover structure characteristics. They also featured a disparate rate of changes in land cover structure.

**Figure 4.** Percentage share of each class.

#### *4.2. Differences in Land Cover*

To achieve a detailed qualitative change analysis, calculations were performed to identify differences in the surface areas according to land cover class for each area between the analyzed years.

Between 1973 and 1951, a drop in the area of land clear of vegetation was observed along the entire analyzed coast fragment with a simultaneous rise in grass-covered and tree-covered terrains. A growth in the water-covered area became noticeable as well. The section of the coast encompassing areas 1–9 constituted an exception to the aforementioned trend. In that terrain, a significant loss of water area of more than 20% was observed, with a simultaneous expansion of the area occupied by other land cover classes (Figure 5).

**Figure 5.** Percentage changes in land cover.

The years 1951–1938 comprised a period during which the tree class cover grew, amounting to 24% in relation to 1938 for area 7 (Figure 5). A simultaneous drop in the area of land covered by grass was recorded. The water-covered area in the individual areas of the entire analyzed coast section fell as well.

The accumulation process in that section was less intensive owing to the currents running alongside the shore and transporting sandy material to the terrains lying further to the west. The lack of beach expansion resulted in limited opportunities for vegetation to inhabit new areas. Here, the dunes were undercut, and there was a narrow strip of grass to the back of the beach, which was followed by a forest. Toward the west, the character of the coast changed as the beach widened, and more initial vegetation was observed, indicating vegetation succession.

As a result of our analyses of the land cover changes occurring between 1938 and 1996, vegetation succession was confirmed for the investigated fragment of the coast, mainly represented by an increasing share of the tree and shrub land cover classes. The first section of the investigated length of coast extended between Swinouj´scie and Mi ˛edzyzdroje (areas ´ 1–15), which had an accumulative character suitable for vegetation succession, and where the most substantial vegetation spread was observed.

The second fragment comprised areas 16 to 23, where slight changes were recorded regarding the surface area occupied by water and an increasing area covered by vegetation. The largest change was observed for the tree class, with an increase of 20% from 1938 to 1996 in area 22.

The third section included areas 24 to 37, where the water surface area had been dominant since 1938, and accounted for approximately 60% of the total land cover. The changes occurring along this fragment of coast involved the expansion of the surface overgrown by trees and shrubs and a decrease in the share of grasses, owing to the lack of opportunities for them to inhabit new grounds. No beach area increments were noted here.

The last fragment of coast comprised areas 38 to 51, where a consistent tendency for the water surface area to expand and a persistent substantial share of the area free of vegetation were observed. A rising share of grass and trees was recorded. This analysis demonstrated that there was a significant correlation between the water surface area and the area free of vegetation, on the one hand, and the changes in the area occupied by vegetation, on the other, which suggested that the structure of the land covered by vegetation was influenced by coastal dynamics.

The relationship between the area covered by water and the area free of vegetation was evident in our subsequent analyses. Where the water surface exhibited rapid changes (10% or more), an inversely proportional change in the area free of vegetation was noticeable.

It is characteristic of accumulative areas that decreases in water-covered terrain and expansion of terrain free of vegetation are accompanied by a growth in the area carrying vegetation. For stabilized shores, tree growth occurs evidently at the expense of grassy areas, with simultaneous small variations in the areas occupied by water and free of vegetation.

#### *4.3. Statistical Analysis of Land Cover Changes*

Finally, the basic fields investigated were analyzed in terms of similarities both in the percentage of the surface occupied by the individual land cover classes and in differences of land cover between the analyzed years. Cluster analysis using Ward's method [46] was implemented, and the results are presented in the tree diagram below (Figure 6).

**Figure 6.** Land cover grouping based on the percentage shares of the individual classes via cluster analysis using Ward's method for all the years in question. Identifier format: YB, where Y—year (A—1938, B—1951, C—1973, D—1996, E—2017), B—basic field number (from 1 to 51).

> Firstly, similarities between the percentage shares of the individual classes were investigated for each year under analysis. According to the measure of the Euclidean

distance, areas closest to one another exhibit the smallest binding distance, which becomes increasingly small as differences between individual areas grow. The above diagram (Figure 6) shows that the areas fell within six characteristic groups. The assignment of the individual areas to the given group in each year was analyzed, and the shares of the particular land cover classes were examined. Furthermore, the average value of the given class in each group was determined. To facilitate the understanding of the spatial distribution of the identified clusters, a diagram was developed that presents the assignment of each single basic field to the specific identified cluster for all the years in question (Figure 7).

**Figure 7.** Diagram showing the assignment of all the areas to groups based on Ward's method.

One characteristic of Group 1 was the approximate 4% and 15% shares of grass and beach surfaces, respectively. Group 2 comprised areas that featured a wide beach, occupying a surface of more than 20%. In Group 3, the approximately 9% share of grass vegetation was a typical occurrence, where it may have been proof of the evolution of a dune that was invaded by grass. Group 4 typically comprised accumulative areas with a slight share of water, a very wide beach, a large proportion of grass in the land cover, and a substantial area taken over by trees. Group 5 was characterized by the approximate shares of 12% and 15% of grass and beaches, respectively, and small areas covered by water. Group 6 was composed of areas where erosion prevailed, grass was virtually non-existent or had a scant share, and the share of trees and shrubs reached approximately 21%.

Based on the spatial distribution of the similarity groups, we were able to identify the variability of the different cluster changes over time. The visible transitions of one type into another in different years constituted a layout that corresponded to the natural coastal structure. In turn, during the analysis of the coast types over time, it was observed that a portion of the areas remained in the same group; for example, area 26 had shown the same coast type since 1938 that was assigned to Group 1. Transitions between the groups could be observed as well. Namely, area 21 once belonged to Group 4 (in 1938), Group 2 (in 1951), Group 1 (in 1973), Group 2 (in 1996), and Group 6 (in 2017), indicating that it experienced significant land cover changes. Overall, 22 areas exhibited the most frequent transformations in structure over time, transitioning between three groups. For 18 areas, a transition between two groups was recorded, with four areas moving between 12 groups.

Furthermore, to better visualize the similarity groups distinguished using Ward's analysis, an average, a minimum, and a maximum of areas were calculated based on all the basic field falling within each cluster. The final averaged land cover characteristic of the clusters is shown in Figure 8, and the extreme values (min. and max.) are shown in Table 2.

**Figure 8.** Average land cover areas in Groups 1–6 classified according to Ward's method.


**Table 2.** Land cover characteristic in Groups 1–6 classified according to Ward's method.

The second part of our cluster analysis, again conducted using Ward's method, focused on similarities in differences of land cover values identified between the individual periods. In the subsequent stage of the study, all the differences observed for the individual periods were subjected to a cluster analysis, and the degree of the relationships between them was examined. The relationships were determined on the basis of the percentage shares of the given classes of land cover in a basic field, much as they were in the case of the classification trees.

Six groups were distinguished based on these results, as shown in Figure 9. The spatial distribution of the land cover share differences between the similarity groups along the coast is shown in Figure 10.

**Figure 9.** Grouping of differences in land cover in the areas based on a cluster analysis using Ward's method for all the areas and all the years in question. Identifier format: YB, where Y—year (A—1973–1996, B—1951–1973, C—1938–1951, D—1996–2017), B—basic field number (from 1 to 51).

**Figure 10.** Areas in Groups 1–6 classified according to Ward's method.

As was the case of the land cover analysis, the similarity classes of land cover changes were described using an average, minimum and maximum of area calculations. An averaged cluster land cover characteristic is shown in Figure 11, and the extreme values (min. and max.) are shown in Table 3.

**Figure 11.** Average differences in land cover areas across Groups 1–6 classified according to Ward's method.

**Table 3.** Land cover differences characteristic in Groups 1–6 classified according to Ward's method.


Finally, the areas assigned to the given group and the structures of the changes in their land cover were verified. In 1938–1973, there were areas assigned to Group 5 where the changes reached a maximum value of 33%. The changes in the structure varied, primarily in terms of the loss of water, accompanied by an increase in areas occupied by beach and vegetation in the form of grass and trees, typical of high accumulation coasts.

In two groups (Groups 2 and 4), the area of land covered by water went up, which was linked to the erosion experienced by those terrains and the accompanying shrinkage of the beach. Group 1 comprised stabilized areas in which a loss of land covered by grass and water was observed, accompanied by a simultaneous incremental increase in the land occupied by trees and sand. Group 2 comprised coastal sections in which a rise in grasscovered land was seen, which may have been evidence of stabilized coastal conditions that benefited the growth of new grass. Increases in grass cover occurred at the expense of beach cover (sand class—areas free of vegetation). Group 3 included very stable coastal areas with low coastal dynamics. Group 4 comprised terrains where accumulation and erosion occurred and was followed by beach expansion and the dwindling of the vegetation cover, as well as terrains with increasing amounts of land covered by water. Group 6 contained territories in which the expansion of the tree cover occurred at the expense of diminishing amounts of land covered by water and beaches. These were the coastal sections along which stable shores experiencing erosion were observed.

The spatial distribution of the results shown in Figure 10 confirms that different areas were assigned to more than one group. According to the analysis, 26 areas were assigned to three different groups, 16 areas to two groups, and 9 to four groups.

The employed statistical analysis demonstrated that the examined fragment of the coast had terrains that were statistically similar to one another despite the distances between them. This finding was also confirmed by the statistical analysis of the differences in land cover changes. The observed area grouping regularity occurred at a distance of approximately 0.5 for both the analyses. A statistical analysis conducted for all the areas in the individual years and for the land cover differences further demonstrated that the grounds had changeable coast characteristics.

The six averaged land cover area groups shown in Figure 8 reflect different types of coast evolution found along the Polish dune coastline, with the spatiotemporal changes of the distinguished types shown in Figures 10 and 11.

#### **5. Discussion**

Aerial photographs are a source of valuable historical information on the vegetation cover and its status [18,19]. As resource managers' and scientists' demands for spatially explicit data continue to grow [18], the use of digitized or digital aerial photographs and the development of automatic analysis techniques can improve the accuracy, consistency, and efficiency of results [47]. Historic information from aerial photographs can also be helpful in monitoring land cover, land cover changes, and coastal dynamics.

Using aerial photos in automated image analysis brings about multiple challenges, including differences in contrast, spectral and spatial resolutions, the solar angle, and the time/year of acquisition. Nevertheless, historical photos in conjunction with plotbased records and appropriate ground truth information provide an important record of vegetation dynamics over time [20].

Investigating land cover and its spatial and temporal changes is essential in coastal sciences. Most recent work deals with their impact on economic development, investment, agricultural policy, and environmental protection policy [33,48–50]. The novelty of the proposed research is that it aims to use land cover changes as a basis for studying coastal dynamics, which allows for the determination of the long-term evolutionary trends of individual fragments of a dune coastline.

This approach differs from many other coastal analyses made worldwide based on simplified global shoreline movement patterns [10,51,52] or site-specific regional analyses of actual coastal changes, that are mostly studied using cross sections of the coast [3,53].

We should keep in mind that, regardless of the line or area indicator methodology, by capturing the historical shoreline position form maps and aerial photographs, one will have to deal with the uncertainty of delineating them. Based on RS data, we are able to mark only the temporal location of a water line. On the one hand, both types of analyses actually neglect the problem of water level change. On the other hand, using area indicators instead of line ones automatically provides additional information about the erosive or accumulative character of change. This appears to be a more reliable method to represent the mean shoreline position for the analyzed period.

Furthermore, the proposed method for determining land cover changes in basic fields containing both water and land surfaces represents an original approach. Although van der Meulen [54] has applied a similar method to the coast of the Netherlands, he only analyzed changes in land cover under the influence of an expected rise in sea level in cells that were 1000-m-long and 50-m-wide, and parallel to the shoreline.

Although the methodology for analyzing coastal dynamics based on basic fields proposed herein has its limitations, mostly associated with the short field, of 1 km in length, the observed changes with regard to the land covered by water, beach, and vegetation do reproduce coastal dynamics with acceptable accuracy and are comparable with earlier research. Thus, the approach where the structure of land cover is analyzed and subsequently used to conduct a statistical analysis of the changes appears to be an efficient technique and an alternative to methods only resorting to cross sections and line indicators [55].

Detailed results for land cover change dynamics confirm the general development tendency demonstrated by the shore as shown by Dudzi ´nska-Nowak [30]. She analyzed the changes of a dune/cliff base occurring between 1938 and 1996 in Rewal Community (research areas no. 32-46, respectively). Based on the analysis, she identified areas subject to morphological erosion and accumulation processes. Her research helped to recognize the variations and dynamics of coast changes. The areas featuring strong accumulation and erosion between 1938 and 1951 turned out either to be stable in the subsequent 1951-73 period or changed their character to either erosive or accumulative. She observed that the period between 1973 and 1996 saw the further disappearance of shore sections

of an accumulative nature in that region. Sections showing an accumulative trend in the previous period transformed into sections of stable shore, whereas sections of stable shore changed their character to erosive [30]. The results of our analysis confirm those tendencies, which proves that the implemented methodology for determining coastal dynamics by analyzing the land cover in basic fields can be compared to line-indicator methods.

Furthermore, the proposed combination of land cover changes derived from historical remote sensing data supported with a statistical analysis is in agreement with the dune coast evolution proposed as part of the general dune shore development model presented by Hesp [56] or Psuty [57]. According to those authors, the coastal foredune is the uppermost and inlandmost component of the sand-sharing system. It has accumulated sand in association with a range of pioneer vegetation types to create a positive landform perched above the dry sand beach. It is the most conservative portion of the coast, undergoing dimensional and temporal changes of a far lesser magnitude and frequency than the sand beach or the offshore zone. In that simple model, the coastal foredune exists on the boundary between the coastal processes on its seaward side and the continental processes occurring landwards. However, many coastal zones are not as simple as this profile, and there are multiple instances of variable dune configurations and areas immediately inland of the dune-beach profile that appear to be morphodynamically related to the processes active in the sand-sharing system [58].

According to our research, performed in the dune areas of the Pomeranian Bay's eastern coast and extending back to 1938, a succession of the vegetation cover was observed, chiefly with regard to the tree class. The fragment of the coast between Swinoujscie and ´ Mi ˛edzyzdroje (areas no. 1 to 15) boasted the most marked vegetation cover succession with simultaneous intensive accumulation processes. Foredunes may be classified into two types: incipient and established. Incipient foredunes are new, or developing, foredunes forming within pioneer plant communities. They may be formed by sand depositing within discrete or relatively discrete clumps of vegetation, or individual plants (types 1a and 1b of Hesp) [56,59]. Such types were identified in areas no. 1 to 15 between 1938 and 1973.

The beaches in the sandy parts of the Polish coast are of considerably varying widths and inclinations, determined by the dynamics in force in the given section of the coast [60]. These findings were confirmed by our results for several of the basic fields analyzed. In areas no. 16 to 23, slight changes to the land area occupied by water were observed, along with substantial proportions of vegetation and beach. The stretch to the east of Dziwnów and to the west of Rewal is dominated by land covered by water (areas 24 to 37). The changes occurring within this fragment of coast are represented by an increasing area carrying tree and shrub land cover classes, which have been successively invading the grass-covered area of land. In turn, a decline in the share of grass cover resulted from its inability to inhabit new terrain, with the beach area failing to expand in that particular section of the coast. In areas 38 to 51, land cover changes demonstrated intensive dynamics. In areas 45 and 50, a substantial growth of grass and tree cover was recorded. In 1938 and 1951, these sections of the coast carried initial vegetation that was characterized by irregular coverage (island coverage). In 1973, a dense structure was observed that covered the dune area, while in 1996 tree succession encroaching onto the area previously occupied by grass was observed. A constant upward trend in terms of land covered by water was accompanied by the persistence of a large proportion of areas free of vegetation, including a wide beach.

Foredunes are very dynamic forms which may rapidly convert into a different form likely to remain in the same location for years to come, until the erosional phase sets in [32,61]. Herein the dynamics of sand class building that forms is an indicator of longterm variability, and could be used as an alternative indicator of coastal dune changes.

The statistical analysis based on Ward's method demonstrated that in terms of land cover structure the studied areas could be classified into six groups according to the shore dynamics presented. Depending on the land cover structure, it is possible to predict the shore evolution stages. These results are consistent with a dune coast evolution model

based on the differences in beach width, a strip of white dunes, and gray dunes overgrown with dense vegetation [62].

#### **6. Conclusions**

The methodology proposed, that accounts for long-term land cover changes as a record of coastal dynamics, is a novel research approach that has never been applied to studies of the Polish Baltic Sea coastal zone. Historical aerial photographs dating back to 1938, 1951, 1973, 1996, and 2017 were processed to conduct a detailed long-term analysis of land cover changes in sections of the dune coast stretching over a total of 51 km. This study was the first for the coastal zone of the Pomeranian Bay's eastern coast, demonstrating that the long-term evolution of dune shores could be analyzed based on changes in land cover using archival RS data.

The methodology developed herein helps discover new possibilities for defining coastal zone dynamics and can be used as an alternative solution to methods only resorting to cross sections and line indicators. The rates of change can be derived from the observed variability of areas of land covered by water or free of vegetation, while the share of the grass and tree classes constitutes a record of the long-term dynamics of the coastal zone. A multidimensional cluster analysis using Ward's method showed that there were six different groups of coast types differing in terms of the land cover structure that is likely to transform over time.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/2072-4 292/13/6/1068/s1, Table S1: Land cover area percentages in each year for areas no. 1–17 [%], Table S2: Land cover area percentages in each year for areas no. 18–36 [%], Table S3: Land cover area percentages in each year for areas no. 37–51 [%].

**Author Contributions:** Conceptualization, field experiments, and formal analysis: A.G.; methodology and data curation: A.G., P.T., P.C., and T.K.; writing original draft, visualization, and validation: A.G. and P.T.; supervision: A.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was co-financed within the framework of the program of the Ministry of Science and Higher Education under the name "Regional Excellence Initiative" during the period 2019–2022; project number 001/RID/2018/19; the amount of financing was PLN 10,684,000.00, and "The Oceanographic Data and Information System eCUDO.pl"—contract no. POPC.02.03.01-00-0062/ 18–00 was funded by the European Union through the Operational Programme Digital Poland 2014–2020 Fund.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy issues.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

