Introducing GEOBIA to Landscape Imageability Assessment: A Multi-Temporal Case Study of the Nature Reserve “K ó zki”, Poland

: Geographic object-based image analysis (GEOBIA) is a primary remote sensing tool utilized in land-cover mapping and change detection. Land-cover patches are the primary data source for landscape metrics and ecological indicator calculations; however, their application to visual landscape character (VLC) indicators was little investigated to date. To bridge the knowledge gap between GEOBIA and VLC, this paper puts forward the theoretical concept of using viewpoint as a landscape imageability indicator into the practice of a multi-temporal land-cover case study and explains how to interpret the indicator. The study extends the application of GEOBIA to visual landscape indicator calculations. In doing so, eight di ﬀ erent remote sensing imageries are the object of GEOBIA, starting from a historical aerial photograph (1957) and CORONA declassiﬁed scene (1965) to contemporary (2018) UAV-delivered imagery. The multi-temporal GEOBIA-delivered land-cover patches are utilized to ﬁnd the minimal isovist set of viewpoints and to calculate three imageability indicators: the number, density, and spacing of viewpoints. The calculated indicator values, viewpoint rank, and spatial arrangements allow us to describe the scale, direction, rate, and reasons for VLC changes over the analyzed 60 years of landscape evolution. We found that the case study nature reserve (“K ó zki”, Poland) landscape imageability transformed from visually impressive openness to imageability due to the impression of several landscape rooms enclosed by forest walls. Our results provide proof that the number, rank, and spatial arrangement of viewpoints constitute landscape imageability measured with the proposed indicators. Discussing the method’s technical limitations, we believe that our ﬁndings contribute to a better understanding of land-cover change impact on visual landscape structure dynamics and further VLC indicator development.


Introduction
The movement of remote sensing software from pixel-based approaches to geographic object-based image analysis (GEOBIA) [1] led to improved workflows for imagery processing, especially land-cover classification [2] and land-cover change detection [3]. GEOBIA has a diverse range of applications [4], including, for example, soil science [5], archaeological research [6], geomorphology [7], forestry [8][9][10], and agriculture [11]. Because of reduced spectral resolution, historical grayscale imagery is usually the object of visual interpretation techniques [12]. Individual [13,14] and sometimes even multiple classes [15] can be digitized on screen by GIS operators. However, these time-consuming tasks can be improved by image segmentation techniques so that the image objects are classified more efficiently. Furthermore, textural analysis enables the supervised classification of single-channel imagery [16] with accuracy comparable to that of on-screen digitalization.
More recently, GEOBIA also became a prerequisite stage for visual landscape indicator development [17,18]. However, the use of GEOBIA for landscape imageability-the ability of a view to make a lasting impression on an observer [19]-and the assessment of visual landscape character (VLC) indicators [19,20] had little investigation to date. Although the contribution of land-use patches to landscape indicator calculation is already a well-recognized topic in landscape ecology [21], landscape quality assessment [22,23], and cultural ecosystem services [24], so far, no comprehensive long-term study exists on the complex relationships between land-cover patch changes and landscape visual structure.
To close this knowledge gap, we put the theoretical concept of using viewpoints as imageability indicators into practice through a multi-temporal land-use case study with three main parts: (i) an investigation of land-cover change over six decades using GEOBIA and high-resolution remote sensing imagery; (ii) deriving imageability indicators using viewpoints and the isovist algorithm; (iii) exploring the relationship between changing land cover and its impacts on the visual landscape. Taking visual landscape [25,26] characteristics as our focus, we interpret the isovist results to determine what they mean for landscape imageability. We hypothesize that the rank and spatial arrangement of viewpoints constitute landscape imageability, while the number and density of viewpoints reveal more about how the structure of a landscape imparts a unique character and even a strong visual image to the landscape user (the observer).

The Importance of Image Segmentation Quality as Prerequisites for Imageability Indicator Calculations
Land-cover patterns, vegetation, and topography are regarded as key landscape features that create a "view" from the perspective of the human observer. Land-cover patterns and landscape features play two main roles in view generation: visual barriers and aesthetic components. We use landscape features as visual barriers to assess VLC, derived through GEOBIA, leaving the more subjective assessment of aesthetic components aside.
Segmentation, as the first stage of GEOBIA [27], is the task of grouping together adjacent pixels that have similar spectral characteristics to define radiometrically homogeneous segments [28]. Each segment represents the core of object-based analysis [1]. As described in Section 3.3, the segments' borders, classified as visual curtains, can be used for imageability indicator calculations.
Although there are several different kinds of segmentation algorithms [29], we chose to use the multi-resolution segmentation (MRS) algorithm [30], which is one of the most widely used and successful [31]. This color-, texture-, shape-, and noise-sensitive algorithm [32] is associated with the difficulty of choosing the optimal scale parameter (SP), as well as the optimal shape and compactness settings. Of the three, however, it is SP which provides the most control of the segmentation process.
The accuracy of automated viewpoint calculations depends on segmentation quality; therefore, finding optimal SP thresholds for each time period of our analysis is essential. However, this can be a challenging task, solved either through a trial-and-error approach [33][34][35] or a more automated method such as using estimation scale parameter (ESP) software [36]. The trial-and-error approach generally omits parameter evaluation and focuses on image objects and their classification for accuracy assessment. However, in the context of our study, where finding an optimal SP for both small and large co-existing landscape features was likely to be challenging and time-consuming, we opted for the more automated approach provided by ESP.
ESP software [37] iteratively generates image objects (segments) at multiple scale levels and calculates the local variance (LV) [38,39] for each of them [40]. The resulting plot of LV change rates (ROC-LV) enables the end-user to select the optimal SP for MRS in the most appropriate manner [41].
The automated scale parameter approach, as an unsupervised method of segmentation accuracy assessment [42][43][44], also provides a more objective basis on which to set SP, a key factor for image object classification accuracy [45,46]. To further improve the process of selecting optimal SP values, an enhancement to ESP software was released as ESP2, an e-Cognition plug-in [36]. The plug-in was successfully used in several studies and is regarded as a credible tool [47][48][49][50] for automating image segmentation at three levels of detail.

Measuring Imageability with the Use of Viewpoints
The VLC indicator theoretical framework, proposed by Ode et al. [19], suggests that "the number and density of viewpoints as imageability indicator(s), could be calculated through visibility analysis using orthophotos or land-cover data and terrain data". Imageability reflects the ability of landscape to create a strong visual image for the observer, thereby making it distinguishable and memorable [19]. Tveit et al. [20] and Ode et al. [51] used imageability to describe a VLC, while Lynch [52] and other urban planners reviewed by Shanken [53] focused on the similar concept of urban legibility.
As reviewed by Ode et al. [19], memorable visual images can be created by iconic elements, landmarks, or land-cover patterns, as well as the number of viewpoints. This theoretical concept does not prejudge whether high or low viewpoint density characterizes imageability. The Ode et al. [19] classification also distinguishes the viewshed size and depth of view as visual scale indicators. Specifically, the visual scale refers to perceptual units of open land delimited by taller vegetation, topography, or human-made objects viewed from a single viewpoint, termed the "landscape room" [20] or "landscape interiors" [54].
Our method uses viewpoints and their properties informed by the VLC imageability indicator classification [19]. We extend the concept of using a viewpoint to measure imageability by distinguishing viewpoint rank and spatial arrangement to describe imageability. To derive viewpoints, we use the isovist algorithm, from which we derived both primary and secondary viewpoints.
First introduced by Benedikt [55], a single isovist is the volume of space visible from a given point, together with a specification of the location of that point (the viewpoint). This behaviorally and perceptually oriented methodology measures what area of space is visible from a predefined viewpoint. One can also think of the isovist as the volume of space illuminated by a point source of light. Every point in physical space has an isovist associated with it. This gives rise to formulating a geometrical model of a minimal isovist set (MIS) covering visible space. In short, the MIS of a simple polygonal region P is the smallest set of viewpoints in P whose union of isovists is equal to P [56]. The two-dimensional MIS is constructed as shown in Figure 1.
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 27 legibility. Therefore, the minimal set of viewpoints and their density are usually used to describe landscape imageability, with their rank and location based on the degree of landscape openness or visual scale. This is achieved with the use of two-dimensional (2D) visibility modelling. gives an overview of most of the landscape and is, therefore, regarded as a primary viewpoint. To see obstructed areas (shaded in red), a secondary viewpoint is necessary (B). However, two small obscured areas remain in red; thus, two more viewpoints (C) are positioned by the MIS algorithm. Full visual coverage of the landscape, therefore, is provided by four viewpoints.

Study Area
Our case study "high natural values" (HNV) grassland landscape belongs to two NATURE 2000 sites: the habitat refuge PLH 14,001 and the bird refuge PLB 140,001 located in the east part of Poland. The central part of this research area (52°21′32″ north (N) 22°51′40″ east (E)) is the Nature Reserve "Kózki" (86.1 ha), located at the Bug river valley, which is also part of Podlasie Bug Gorge Landscape park. This relatively flat area is located at mean 117 m above sea level, 2 km from the town of Siemiatycze (Figure 2). Its high natural values (HNVs) are created by expanses of gray and white dunes overgrown by xeric sand calcareous grassland communities [62]. Further grassland gives an overview of most of the landscape and is, therefore, regarded as a primary viewpoint. To see obstructed areas (shaded in red), a secondary viewpoint is necessary (B). However, two small obscured areas remain in red; thus, two more viewpoints (C) are positioned by the MIS algorithm. Full visual coverage of the landscape, therefore, is provided by four viewpoints.
The MIS, where viewpoints have a high field of view, are termed primary viewpoints, whereas secondary viewpoints have relatively low fields of view and are usually more numerous. In the context Remote Sens. 2020, 12, 2792 4 of 25 of VLC, specifically imageability, the primary viewpoints create a sense of the infinite and mystery of what lies in the landscape beyond the vantage point [57]. In general, long-distance views relate to our ancestors' experiences of open landscapes, for which the possibility of looking at a distance determines the safety and survival in the living environment. On this basis, a long-distance view provided by the primary viewpoints contributes to landscape users' psychological comfort, sense of personal freedom, security, and mental pleasure [58], as well as contemplative experience of the space [59]. According to the most recent contemplative landscape research findings that used neuroscience tools to evaluate the health impact of landscape exposure [60], the depth of the view and the visibility of landscape layers (fore-, middle-, and background) is one of the key factors contributing to the restorative effect in the human brain [61]. Secondary viewpoints, meanwhile, tend to entice an observer toward deeper fields of view. Their contribution to imageability is less important unless they form a characteristic system (e.g., maze), constituting a specific landscape legibility. Therefore, the minimal set of viewpoints and their density are usually used to describe landscape imageability, with their rank and location based on the degree of landscape openness or visual scale. This is achieved with the use of two-dimensional (2D) visibility modelling.

Study Area
Our case study "high natural values" (HNV) grassland landscape belongs to two NATURE 2000 sites: the habitat refuge PLH 14,001 and the bird refuge PLB 140,001 located in the east part of Poland. The central part of this research area (52 • 21 32" north (N) 22 • 51 40" east (E)) is the Nature Reserve "Kózki" (86.1 ha), located at the Bug river valley, which is also part of Podlasie Bug Gorge Landscape park. This relatively flat area is located at mean 117 m above sea level, 2 km from the town of Siemiatycze (Figure 2). Its high natural values (HNVs) are created by expanses of gray and white dunes overgrown by xeric sand calcareous grassland communities [62]. Further grassland community details of "Kózki" Nature Reserve were described by Warda et al. [63] and Kulik et al. [64].
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 27 legibility. Therefore, the minimal set of viewpoints and their density are usually used to describe landscape imageability, with their rank and location based on the degree of landscape openness or visual scale. This is achieved with the use of two-dimensional (2D) visibility modelling. gives an overview of most of the landscape and is, therefore, regarded as a primary viewpoint. To see obstructed areas (shaded in red), a secondary viewpoint is necessary (B). However, two small obscured areas remain in red; thus, two more viewpoints (C) are positioned by the MIS algorithm. Full visual coverage of the landscape, therefore, is provided by four viewpoints.

Study Area
Our case study "high natural values" (HNV) grassland landscape belongs to two NATURE 2000 sites: the habitat refuge PLH 14,001 and the bird refuge PLB 140,001 located in the east part of Poland. The central part of this research area (52°21′32″ north (N) 22°51′40″ east (E)) is the Nature Reserve "Kózki" (86.1 ha), located at the Bug river valley, which is also part of Podlasie Bug Gorge Landscape park. This relatively flat area is located at mean 117 m above sea level, 2 km from the town of Siemiatycze ( Figure 2). Its high natural values (HNVs) are created by expanses of gray and white dunes overgrown by xeric sand calcareous grassland communities [62]. Further grassland community details of "Kózki" Nature Reserve were described by Warda et al. [63] and Kulik et al. [64].

Materials and Methods
Our methodology links GEOBIA with VLC assessment procedures. The GEOBIA stages include image preprocessing, segmentation, classification, and accuracy assessment. The GEOBIA

Materials and Methods
Our methodology links GEOBIA with VLC assessment procedures. The GEOBIA stages include image preprocessing, segmentation, classification, and accuracy assessment. The GEOBIA classification result is then input as a two-dimensional landscape model for viewpoint calculation using the MIS algorithm implemented in the Isovists (version 2.3) [65] software. We use the derived viewpoints to estimate the minimum set and measure their density, spacing, and hierarchy-all core indicators of landscape imageability [19].

Remote Sensing Imagery Pre-Preprocessing
The relatively small, low-contrast land-cover patches of our case study site make the use of VHR images-the most suitable for the case-study purposes. For this reason, global Earth Observation (EO) programs with coarser pixels (like Landsat NASA or Sentinel-2 ESA) were excluded in favor of historical aerial, VHR satellite, and UAV imagery. An exception was CORONA KH-4 declassified imagery, which, despite its relatively coarse spatial resolution, allowed us to include a time period from the 1960s. Historical imagery was initially sourced from the Polish National Geodetic and Cartography Repository (PNGCR). Only imagery from the growing season and without cloud cover was considered. Off-nadir imagery was included to avoid the exclusion of potentially useful imagery.
To create a more complete time series, we also obtained more recent imagery from satellite (Pleiades-1B) and UAV platforms. The resulting collection of imagery had varied characteristics, including differing scale or ground sample distance (GSD) and spectral resolution (from grayscale to four-band compositions). The complete list of remote sensing datasets used in this study, along with basic technical characteristics, is listed in Table 1 and presented in Figure 3.

Remote Sensing Imagery Pre-Preprocessing
The relatively small, low-contrast land-cover patches of our case study site make the use of VHR images-the most suitable for the case-study purposes. For this reason, global Earth Observation (EO) programs with coarser pixels (like Landsat NASA or Sentinel-2 ESA) were excluded in favor of historical aerial, VHR satellite, and UAV imagery. An exception was CORONA KH-4 declassified imagery, which, despite its relatively coarse spatial resolution, allowed us to include a time period from the 1960s. Historical imagery was initially sourced from the Polish National Geodetic and Cartography Repository (PNGCR). Only imagery from the growing season and without cloud cover was considered. Off-nadir imagery was included to avoid the exclusion of potentially useful imagery.
To create a more complete time series, we also obtained more recent imagery from satellite (Pleiades-1B) and UAV platforms. The resulting collection of imagery had varied characteristics, including differing scale or ground sample distance (GSD) and spectral resolution (from grayscale to four-band compositions). The complete list of remote sensing datasets used in this study, along with basic technical characteristics, is listed in Table 1 and presented in Figure 3.

Imagery Processing Method
Grayscale images (1957,1965,1973) were retrieved as raw 8-bit raster files that needed to be georeferenced and mosaicked into a single image. Similar to other cases of historical imagery pre-processing [66], metadata like interior and exterior camera orientation (applied to aerial imagery) were unknown. Therefore, we applied the image to an image co-registration technique [14,[67][68][69], as an alternate method known to yield satisfactory results especially in low-relief terrains [70]. A highresolution (GSD 0.25 m, horizontal accuracy 0.15 m) contemporary orthoimage [71] was used as a reference to conduct co-registration.
For each georeferenced image, at least 10 equally distributed ground control points (GCPs) were manually selected through careful inspection of the reference and archival images. The task of identifying the ground features was conducted by an analyst with good knowledge of the study area from field surveys. As GCPs, we used identifiable, stable ground features such as road intersections, buildings, bridges, and eventually natural features (e.g., single trees), which can also be used for GCP in the absence of anthropogenic features [72]. The height (z-value) of each GCP was read from the SRTM elevation model (DEM Version 4) to enhance the polynomial transformation of the interpolator [73,74]. Furthermore, because we used imagery spanning 61 years of landscape change, we selected imagery-specific GCPs for each period.
Image co-registration was done to single-pixel accuracy. A third-order polynomial transformation that resulted in the lowest RMSE and least visible image distortion was selected for image co-registration (reported in Table 1 as metadata). All datasets were registered in the EPSG 2180 coordinate system. Co-registered images were then mosaicked into single raster datasets for each time period and clipped to the extent of the 2018 UAV imagery. The UAV imagery encompassed the extent of the "Kózki" nature reserve, but was otherwise the imagery with the smallest spatial extent of our study ( Figure 3H); thus, it became our study area boundary. Finally, all archival imagery (except CORONA) was down-sampled to the same GSD (0.5 m) to make the comparison among the imagery easier in subsequent processing steps.
For the imagery retrieved from the PNGCR as orthoimages, no pre-processing was done except for clipping to the research area boundary and resampling to 0.5 m. The Pleiades orthoimages were delivered from Airbus Space as two (pan, multispectral) Tiff files after atmospheric correction (scene ID: DS_PHR1B_201507040930594_FR1_PX_E022N52_1008_03164). To achieve the GSD adopted in our study, the Pleiades imagery was pansharpened using the Gram-Schmidt method [75] at band weights dedicated for this particular sensor (R 0.166, B, G 0.167, IR 0.5) [76]. The most recent orthophoto was acquired on 26 June 2018 from a UAV platform provided by the University of Wroclaw (Department of Geoinformatics and Cartography; Poland). According to the UAV service provider report [77], the flight was done with the use of an eBee drone equipped with a Canon S100 camera at an ATO height of 148 m AGL. Orthoimagery was produced by the UAV contractor using Metashape (Agisoft). In total, 36 signalized and 23 natural GCPs were used, with horizontal and vertical accuracies of 0.05 m. Finally, the delivered RGB and CIR imagery was merged into four-band Tiff files and resampled from GSD 0.39 to 0.5 m, as with all other source data.

Segmentation and Segment Evaluation Method
The multi-resolution segmentation (MRS) algorithm [30] implemented in e-Cognition Developer (Trimble Geospatial) software was applied for imagery segmentation. For shape and compactness segmentation parameters, the default values of 0.1 and 0.5, respectively, were retained. However, for the key factor of segmentation, the spatial parameter (SP) from an ROC-LV analysis was used.
The ROC-LV graph was generated as a bottom-up process of 150 loops with the use of ESP2 software. For each set of imagery, the ROC-LV was analyzed to select five SP candidates. Using this unsupervised method helped to avoid inefficient, step-by-step SP candidate testing [78]. To decide which among the five SP candidates was the most appropriate, a supervised method of assessing segmentation accuracy was applied [79]. Lucieer and Stein [80,81] pioneered the method of using topological object metrics to assess under-and over-segmentation; thus, segmentation accuracy metrics are now well documented [33,79,82]. Specifically, GEOBIA segments (GSs) are compared to ground-truth reference segments (RSs). Segmentation accuracy metrics can be calculated when several RSs match with one or more GS in terms of shape and size [78]. RSs are usually digitized with the use of orthoimagery [83]. In GEOBIA, the well-defined shape of arable fields [33,84], water bodies [85], trees [36], and some anthropogenic structures (e.g., building footprints) are preferred as RSs.
All RSs were digitized manually (on-screen vectorization) by two independent GIS operators. Next, in accordance with Zhang et al. [86], each RS was discussed by the GIS operators, and only those agreed upon were used for segmentation accuracy metric calculation. We also digitized each time period's RS segments separately. By doing this, the land-cover changes, vegetation phenology. and spectral properties of each imagery dataset were captured by each RS set. Furthermore, as recommended by Moller et al. [33] and Marup et al. [87], the RS count and area class proportions [79] were also restricted. We used a set of 50 RSs at area class proportions set to 3:1:1 for small (≤0.09 ha), medium (0.1-0.2 ha), and large (>0.2 ha) RS area classes, respectively. An exception was made for the 1965 CORONA imagery, where its lower spatial resolution limited the number of small RSs we were able to confidently digitize.
We chose a higher proportion of small RSs due to the fact that the optimal SP was expected to create segments small enough to catch homogeneous patches of shrubs growing over the open sands of the study area, but that also resulted in over-segmentation of the neighboring grassland patches. Because of strong shadow effects (especially in 2010 and 2018), mixed forest areas could not be represented as a single segment, with pre-testing showing that homogeneous segmentation of these areas was difficult even at SPs greater than 150, as well as at ESP2 Level 3.
To assess segmentation accuracy, we used four common metrics: (1) f-score (FS) metrics combined from precision (p) and recall (r) measures [42], (2) the segmentation error (SE) [88], (3) the Jaccard index (JI) [89,90], and (4) accuracy (AC). The applied metrics' mathematical formulae are presented in Equations (1)-(4). The metrics range in value from 0-1, with lower SE values showing higher RS to GS similarity and, thus, the optimal SP, with the inverse being true for FS, JI, and AC. To select optimal values, SPs with the best metric scores were selected for each imagery date. where (recall), n is the number of GS segments, RS is the reference segment, GS i is the i-th GEOBIA segment and E i = envelope(RS ∪ GS i ).

Land-Cover Class Nomenclature
For our multi-temporal analysis, five land-cover classes were used. According to Anderson's [91] land-cover classification levels (LCCL) used by the USGS, panchromatic images are only suitable for level I (LCCL-1) classification because of their low spectral resolution. However, from the perspective of landscape visual property analysis, only two of the nine LCCL-1 classes (built-up and forest) can be used as visual curtains. On the other hand, the sub-meter spatial resolution of archival RS imagery allows for the extraction of other important land-cover classes at LCCL-2 (e.g., shrubs, orchards, vineyards, and bare exposed rock). Furthermore, Corine Land Cover (CLC) nomenclature provides additional class description details that can affect imageability indicator calculation. For imageability analysis, the appointed LCCL should include compulsory land-cover classes (c-classes) which create visual curtains and, optionally, case-study specific classes (s-classes), including, in this case, xeric sandy grassland, open sand, and water. The adopted classification nomenclature is presented in Table 2, along with references to USGS [91] and CLC nomenclature. Table 2. Land cover classes used in the study. Descriptions are provided, along with references to USGS and Corine Land Cover (CLC) nomenclature. Only references that most accurately describe the case-study land cover are listed (e.g., CLC Level 1 forest is not reported because, at Level 1, forest is merged with semi-natural areas). LCCL-land-cover classification level.

GEOBIA Classification Methodology
Taking advantage of image object geometry, texture, and brightness, we adopted the Reference [67] framework of combining manual and supervised image object classification, along with post-processing procedures. Firstly, water bodies and building footprints were manually digitized (on-screen visual interpretation) and incorporated as classified segments into an MRS workflow. This helped to solve the problems of water body misclassification caused by reflection and the presence of sandbanks, as well as misclassification of buildings with gabled roofs due to high contrast among roof surfaces. All other segments generated at optimal SP (Section 3.3) were classified by supervised methods incorporating training (ts) and validation (vs) ground-truth samples.
For each imagery period, a set of 250 land-cover samples was created by a GIS operator, with 50 in each land-cover class. Only well-recognized samples were selected with the use of time series orthoimages, as well as knowledge from previous field work and information on the vegetative community structure of the case study area [62]. For this reason, a random sampling method for ground-truth sample selection was not applied. Ground-truth samples were divided into two equal sets of ts used for classifier training and vs used for accuracy assessment.
The classification procedure was conducted with the use of e-Cognition Developer as an automated process, along with a classification error matrix [92] and accuracy assessment summary including overall, user, and producer accuracy and Kappa measure [93]. Automation of the classification process allowed us to test several of the most popular classifiers implemented in e-Cognition: random forest (RF) [94], support vector machine (SVM) [95,96], and K-nearest neighbor (KNN) [97]. Each classifier received the same input parameters: brightness, degree of skeleton branching, curvature, and textures after Haralick [98] (GLCM homogeneity, entropy, mean, dissimilarity). The best classifier was chosen based on confusion matrix results [92]. Classification post-processing included identification of misclassified image objects and manual correction. The main focus of post-classification correction was for c-class objects, as they were the most important for imageability analysis. Finally, the c-class GEOBIA results were converted to CAD format (.dxf) to create land-cover edges (LCEs) that became visual barriers for MIS analysis (Section 3.5).

Isovist and Imageability Indicator Method
To assess the visual consequences of landscape change, three imageability indicators were calculated: number of viewpoints (Vn), viewpoint density (VD), and viewpoint spacing (VS) (Equations (5) and (6)). The role of these metrics is to objectively describe VLC changes over the 61 years of our study, rather than subjectively landscape aesthetics. Therefore, the meaning of the imageability indicator output is descriptive, rather than normative, and it is not meant to comment on landscape quality. In general, low indicator values characterize landscapes with a strong visual impression created by openness. Gradual increases in VLC metric values indicate a shift to a more closed and less visually connected landscape, with visual impressions resulting more from landscape rooms (interiors) rather than a deep field of view. VD = Vn/A, where Vn is the number of viewpoints, and A is the research area.
The indicator values rely on the viewpoint number (Vn) calculated during MIS analysis implemented in ISOVISTS software [65]. To initiate the MIS analysis, a strategic point was placed common to the LCE generated for each time period of our analysis. In our case, this was the middle of Bug River bridge in the northeastern part of the study area. Next, a series of minimally overlapping isovist fields was generated until the entire analysis area was covered. For each isovist, the points of maximum permitted field and minimum occlusivity are marked as viewpoints, connected by at least one uninterrupted line of sight. The algorithm marked the viewpoints of the dominant isovist [99] as large rings proportional to the area of each viewpoint's isovist field. To better describe the imageability of each landscape, we then classified viewpoints into primary (the most prominent viewpoints) and secondary viewpoints, and we cartographically represented their fields of view. Additional conclusions characterizing the changing VLC were formulated based on the spatial arrangement of viewpoints.

Results
The study returned three types of multi-temporal results: land-cover maps, viewpoint maps, and imageability indicators. The accuracy of the land-cover segmentation and classification is reported in Sections 4.1-4.3. Section 4.4 presents the results of the imageability analysis, including the number, density, spatial distribution, and hierarchy of viewpoints regarded as imageability indicators. This is followed by our discussion of how changes in land cover affect landscape imageability, and how imageability can be described using viewpoint properties as core indicators.

SP Candidate Results
Scale signatures are presented as ROC-LV graphs (Figure 3) with scale ranges from 1-150. Over the study periods, the maximum value of LV ranged from 15.05 (1973) to 72.43 (2015), with an average across most years closer to 20. This indicates a degree of similarity for object heterogeneity, as well as a similar spatial dimension of the most characteristic objects across time periods. The LV values in 2015 are outliers to this trend, with high LV values possibly the result of the pan-sharpening procedure used in that year. All LV graphs show a rising trend, which confirms the similarity of the compared areas.
Because of the smoothed shape of the LV line, ROC lines were used to pick SP candidates. Local variation in the ROC line made simple selection of thresholds difficult; thus, the vertical ROC axes were rescaled to a value of 2 to enable better interpretation of prominent peaks in the ROC line against the x-axis. The most significant scale levels were identified for the grayscale CORONA imagery, which provides a good example for our selection method. While an SP value of 107 was identified as optimate by the automatic ESP2 analysis at Level 1, a more significant threshold was identifiable in the ROC-LV graph at SP 120, which was our final selection for segmentation (Figure 4). The following sets of SP were selected for segmentation using MRS in each year: 45, 75, 91, 135, 149 (1957)  The ROC-LV graph analysis resulted in 40 SP candidates with scale ranges from 23 to 149 to be used as MRS input parameters. The comparison of subsequent segmentation results revealed that the relationships between SP values and segment properties differed depending on the spectral properties and quality of imagery used. The same SP value applied in segmentation of two different images of the same research area results in a different image object count, shape, and size. In our case study, five of the 40 SPs were repeated for at least two segmentation procedures. On the one hand, repeated SPs reveal some common geoprocessing properties (e.g., comparable minimum segment size); however, on the other hand, different segmentation results are due to differences in image spectral characteristics and quality, and, above all, in land-cover changes between compared landscapes. The summary statistics presented in Appendix A (Table A1) confirm the need for individual selection of SPs for each imagery dataset, even for the same research area.

The Reference Segment Digitalization results
This also implies that MRS result evaluation should be carried out with the use of RS and shape The ROC-LV graph analysis resulted in 40 SP candidates with scale ranges from 23 to 149 to be used as MRS input parameters. The comparison of subsequent segmentation results revealed that the relationships between SP values and segment properties differed depending on the spectral properties and quality of imagery used. The same SP value applied in segmentation of two different images of the same research area results in a different image object count, shape, and size. In our case study, five of the 40 SPs were repeated for at least two segmentation procedures. On the one hand, repeated SPs reveal some common geoprocessing properties (e.g., comparable minimum segment size); however, on the other hand, different segmentation results are due to differences in image spectral characteristics and quality, and, above all, in land-cover changes between compared landscapes. The summary statistics presented in Appendix A (Table A1) confirm the need for individual selection of SPs for each imagery dataset, even for the same research area.

The Reference Segment Digitalization Results
This also implies that MRS result evaluation should be carried out with the use of RS and shape similarity metrics. Examples of RS polygon digitization are presented in Figure 5B-D, along with overlapping RS polygon results (A). The overlap count number identifies zones of clearly distinguishable homogeneous patches. In the southern part of the study area, they are arable field shapes, and, in the central part, they are water bodies (or patches of aquatic vegetation), while, in the northcentral section, RS shapes correspond more to open sand or groups of trees, characteristic of the case study area. The mean size of each set of 50 RS was 0.13, 0.41, 0.12, 0.11, 0.10, 0.13, 0.12, and 0.09 for 1957, 1965, 1973, 1997, 2006, 2010, 2015, and 2018, respectively. The 1965 outlier is a result of low CORONA imagery resolution, which was also the reason for the methodological choice of using equal proportions of small, medium, and large RSs for this year.
Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 27 low CORONA imagery resolution, which was also the reason for the methodological choice of using equal proportions of small, medium, and large RSs for this year.

Segmentation Accuracy Results
SP is the main parameter that controls MRS and determines the size and shape of resulting segments. To know which SP value provides the greatest similarity between GSs and RSs, we use shape similarity metrics, as well as precision and recall, to compare segmentation quality. The precision and recall curves calculated for five SP candidates are shown in Figure 6, while extensive details on segmentation goodness metrics are provided in Appendix A (Table A2). Each curve of Figure 6 has a falling trend from upper left to the lower right, typical of these opposing indicators, with curves located in the upper-right part of the graph indicative of high segmentation quality [86]. The worst segmentation performance is for the CORONA imagery (1965); however, due to that imagery's low GSD, these results are not unexpected. At the same time, overall image quality is important, with the 1957 aerial imagery achieving the worst segmentation quality. The curves of 2006 (green) and 2010 (light blue), and, to a lesser extent, the curves of 2018 (purple), 1997 (light green), and 1973 (yellow) show higher segmentation quality for the first three SPs. The central location of the 2015 curve (blue) highlights the differences between MRS performed on satellite

Segmentation Accuracy Results
SP is the main parameter that controls MRS and determines the size and shape of resulting segments. To know which SP value provides the greatest similarity between GSs and RSs, we use shape similarity metrics, as well as precision and recall, to compare segmentation quality. The precision and recall curves calculated for five SP candidates are shown in Figure 6, while extensive details on segmentation goodness metrics are provided in Appendix A (Table A2). Each curve of Figure 6 has a falling trend from upper left to the lower right, typical of these opposing indicators, with curves located in the upper-right part of the graph indicative of high segmentation quality [86]. The worst segmentation performance is for the CORONA imagery (1965); however, due to that imagery's low GSD, these results are not unexpected. At the same time, overall image quality is important, with the 1957 aerial imagery achieving the worst segmentation quality. The curves of 2006 (green) and 2010 (light blue), and, to a lesser extent, the curves of 2018 (purple), 1997 (light green), and 1973 (yellow) show higher segmentation quality for the first three SPs. The central location of the 2015 curve (blue) highlights the differences between MRS performed on satellite imagery compared with lower-altitude imagery.

Segment Classification Accuracy Results
The MRS calculated for each of the eight analyzed time-frame imageries resulted in individual image object segments, which were the object of further supervised classification (excluding water and building footprints processed manually). We were interested in merging all segments belonging to particular land-cover classes and then extracting LCE as input for imageability indicator calculation. The results of imageability analysis depend on object classification accuracy; since the research is focused more on imageability rather than GEOBIA classification efficiency, the details on contributions of the different variables to classification accuracy are limited to accuracy summary statistics. Using 125 training samples (25 per class), various levels of segment classification accuracy were obtained depending on the classifier used. The accuracy measures (overall accuracy, as well as user and producer accuracy) as a result of 24 confusion matrix calculations, as well as additional Kappa statistics, are presented in Appendix B (Table A3). The resulting overall accuracy (OA) ranged from 92% to 46%. In general, the results obtained at OA above 85% were accepted as sufficient for further imageability analysis. The highest accuracy was obtained in case of Pleiades image classification, where the RF classification resulted in OA at 92%. Furthermore, the RF proved to be an effective algorithm in most classified segments. In two of 24 classifications, KNN (1965) and SVM (1973) resulted in higher OA than RF (OA 73% and 85% for 1965 and 1973, respectively). Furthermore, both c-classes were classified at a producer accuracy not worse than 88% for forest segments and 76% for shrub segments. Moreover, segments of sand (an r-class) as the brightest part of imagery were correctly classified at a producer accuracy between 100% and 84%, except for 1997, where the producer accuracy was 60%. The Xeric sand grassland was the most often misclassified Although informative, changes in segmentation quality rates reflected by precision and recall curves cannot indicate optimal SP values. Therefore, additional segmentation quality metrics were used (FS, SE, JI, and AC), the results of which are reported in detail in Appendix A (Table A2), along with the selected optimal SP values. In the case of 1957, an SP of 75 was chosen due to its highest values of FS, JI, and AC, as well as the lowest SE value. Furthermore, the SP of 75 turned out to be the only SP value also confirmed by the automated ESP2 results as optimal for Level 1 segmentation (Figure 4a). In other cases, the choice of an optimal SP was confirmed by the compliance of all (2006, 2010, and 2018) or at least three of four (1965, 1973, and 1997) calculated metrics (details in Table A2, Appendix A). Finally, the following SPs were selected as optimal for the case study MRS: 75, 72, 72, 62, 88, 108, 97, and 92 for 1957, 1965, 1973, 1997, 2006, 2010, 2015, and 2018, respectively.

Segment Classification Accuracy Results
The MRS calculated for each of the eight analyzed time-frame imageries resulted in individual image object segments, which were the object of further supervised classification (excluding water and building footprints processed manually). We were interested in merging all segments belonging to particular land-cover classes and then extracting LCE as input for imageability indicator calculation. The results of imageability analysis depend on object classification accuracy; since the research is focused more on imageability rather than GEOBIA classification efficiency, the details on contributions of the different variables to classification accuracy are limited to accuracy summary statistics. Using 125 training samples (25 per class), various levels of segment classification accuracy were obtained depending on the classifier used. The accuracy measures (overall accuracy, as well as user and producer accuracy) as a result of 24 confusion matrix calculations, as well as additional Kappa statistics, are presented in Appendix B (Table A3). The resulting overall accuracy (OA) ranged from 92% to 46%. In general, the results obtained at OA above 85% were accepted as sufficient for further imageability analysis. The highest accuracy was obtained in case of Pleiades image classification, where the RF classification resulted in OA at 92%. Furthermore, the RF proved to be an effective algorithm in most classified segments. In two of 24 classifications, KNN (1965) and SVM (1973) resulted in higher OA than RF (OA 73% and 85% for 1965 and 1973, respectively). Furthermore, both c-classes were classified at a producer accuracy not worse than 88% for forest segments and 76% for shrub segments. Moreover, segments of sand (an r-class) as the brightest part of imagery were correctly classified at a producer accuracy between 100% and 84%, except for 1997, where the producer accuracy was 60%. The Xeric sand grassland was the most often misclassified community, usually confused with agricultural background, especially in the case of SVM. The producer accuracy of the best Xeric sand grassland classification ranged from 40% (1965) to 100% (2015). Despite this, the Xeric sand grassland, as an r-class, does not impact further visibility analysis. From the perspective of the case study, the Xeric sand grassland patch locations relative to viewpoints have key importance for analyzing HNV. Finally, the misclassified segments were the object of manual post-classification to provide reliable land-cover results. The final land-cover maps, along with segmentation results, are presented at Figure 7.

Land Cover Changes
Main changes within the r-class referred to open sand patches; thus, characteristics for this open landscape decreased by 79.6% (from 20.23 ha to 4.13 ha). The results of land-use change dynamics presented as a percentage of land-cover class area revealed the increase in forest area from 8.97 ha (1957) to 80.07 ha (2018) and relatively stable area of the shrub community ranging from 3.44-7.70 ha, and most recently 3.89 ha in 2018 (Figure 8). Only noticeable increases in the shrub community area between 1997 and 2006 are the result of natural succession, which, in subsequent years, was eliminated under the active sheep grazing program for the "Kózki" nature reserve in 2010 [64]. This trend is also reflected in patch number (NP) and density (PD) analysis results (Table 3). Taking into account only c-class patches as the most relevant for VLC, the landscape pattern changed from small to large patches, maintaining their NP at a relatively unchanged level (NP oscillation between 465 and 598, excluding coarse 1965 imagery). These land-cover change directions reflect a question regarding VLC transformation. To answer this, the results on imageability indicators are presented in the next section.    [64]. This trend is also reflected in patch number (NP) and density (PD) analysis results (Table 3). Taking into account only c-class patches as the most relevant for VLC, the landscape pattern changed from small to large patches, maintaining their NP at a relatively unchanged level (NP oscillation between 465 and 598, excluding coarse 1965 imagery). These land-cover change directions reflect a question regarding VLC transformation. To answer this, the results on imageability indicators are presented in the next section.

The Imageability of Changing Landscape Interpretation
The MIS calculations resulted in unique viewpoint arrangements ( Figure 9A-H) and indicator values ( Table 4) that reflect the VLC changes. Specifically, 1957, as a starting point of 61 years of transformation, displayed its openness, limited only by the distant forest patches. Three primary visibility points created a 1.32-km route, which allows admiring the overviewed landscape ( Figure  9A), whose character is created by open sands and unobstructed water bodies. These characteristics are also confirmed by the lowest imageability indicator values, VD = 62 point/km 2 and VS = 0.13 km. In the case of 1965, due to no shrubs detected as view curtains, only two primary viewpoints provide the central section of the research area overview. Over time, this unique visual character changed irreversibly. In 1973, due to afforestation, the new primary viewpoint was created and maintained until 2018 within the Bug River northern section. Two more primary viewpoints were delimited in the central zone; however, their location favors close-range observation of open sands and Xeric sand grassland, constituting the HNV of the analysis. In 1997, the value of VS became less than 0.1 (0.09), which, together with viewpoint location analysis, can be interpreted as imageability created more by the sequences of landscape interiors and view panoramas [54] rather than stunning openness. Importantly, 1997 and 2006 gradually closed the visual connectivity with the central zone of the research area. This visual connectivity was maintained until 2010, whereas, in 2015, the central section became visually separated (Figure 9 G,H). Visual connectivity closing was accomplished by an increase in viewpoint number with gradual predominance of secondary viewpoints. Furthermore, as the number of viewpoints rises, their field of view is reduced, which underlines the role of close-range visual impressions in imageability creation. Patches of open sands and Xeric sand grassland, as well as the Bug River natural valley, constitute the HNV of the contemporary "Kózki" reserve are. From 2006 until 2018, grassland communities were the location of primary viewpoints. Importantly, the viewpoint sequences created a visibility route from the east (the road) to the west, guiding the landscape user through sand and grassland landscape to border river. The imageability of the contemporary reserve "Kózki" is created by a sequence of landscape interiors closed by forest walls and sliced by spared scrub communities. Concluding, the interpretation of the indicator values allows capturing a very general VLC (open, closed, visually connected/disconnected), which provides the basis for changing landscape imageability description only in combination with viewpoint spatial distribution.

The Imageability of Changing Landscape Interpretation
The MIS calculations resulted in unique viewpoint arrangements ( Figure 9A-H) and indicator values ( Table 4) that reflect the VLC changes. Specifically, 1957, as a starting point of 61 years of transformation, displayed its openness, limited only by the distant forest patches. Three primary visibility points created a 1.32-km route, which allows admiring the overviewed landscape ( Figure 9A), whose character is created by open sands and unobstructed water bodies. These characteristics are also confirmed by the lowest imageability indicator values, VD = 62 point/km 2 and VS = 0.13 km. In the case of 1965, due to no shrubs detected as view curtains, only two primary viewpoints provide the central section of the research area overview. Over time, this unique visual character changed irreversibly. In 1973, due to afforestation, the new primary viewpoint was created and maintained until 2018 within the Bug River northern section. Two more primary viewpoints were delimited in the central zone; however, their location favors close-range observation of open sands and Xeric sand grassland, constituting the HNV of the analysis. In 1997, the value of VS became less than 0.1 (0.09), which, together with viewpoint location analysis, can be interpreted as imageability created more by the sequences of landscape interiors and view panoramas [54] rather than stunning openness. Importantly, 1997 and 2006 gradually closed the visual connectivity with the central zone of the research area. This visual connectivity was maintained until 2010, whereas, in 2015, the central section became visually separated ( Figure 9G,H). Visual connectivity closing was accomplished by an increase in viewpoint number with gradual predominance of secondary viewpoints. Furthermore, as the number of viewpoints rises, their field of view is reduced, which underlines the role of close-range visual impressions in imageability creation. Patches of open sands and Xeric sand grassland, as well as the Bug River natural valley, constitute the HNV of the contemporary "Kózki" reserve are. From 2006 until 2018, grassland communities were the location of primary viewpoints. Importantly, the viewpoint sequences created a visibility route from the east (the road) to the west, guiding the landscape user through sand and grassland landscape to border river. The imageability of the contemporary reserve "Kózki" is created by a sequence of landscape interiors closed by forest walls and sliced by spared scrub communities. Concluding, the interpretation of the indicator values allows capturing a very general VLC (open, closed, visually connected/disconnected), which provides the basis for changing landscape imageability description only in combination with viewpoint spatial distribution.

Discussion
This study extends the application of GEOBIA by utilizing landscape imageability indicators to assess VLC, both qualitatively (VLC description) and quantitatively (indicator values). We show that GEOBIA results-specifically, the use of classified objects (forest, shrubs, and building edges) as visibility curtains-can support VLC assessment, along with additional landscape imageability indicators. To the best of our knowledge, the idea of using viewpoints as imageability indicators was, to date, only presented as a theoretical concept [19]. Therefore, a discussion of indicator value results cannot be conducted with respect to other studies. Our segmentation accuracy results, meanwhile, do correspond with similar studies [100,101], although it is often a challenge to compare segmentation accuracies among GEOBIA case studies due to unique methods, remote sensing imagery, and landscape types.
Although there are no studies to compare this one to directly, we still identified some important

Discussion
This study extends the application of GEOBIA by utilizing landscape imageability indicators to assess VLC, both qualitatively (VLC description) and quantitatively (indicator values). We show that GEOBIA results-specifically, the use of classified objects (forest, shrubs, and building edges) as visibility curtains-can support VLC assessment, along with additional landscape imageability indicators. To the best of our knowledge, the idea of using viewpoints as imageability indicators was, to date, only presented as a theoretical concept [19]. Therefore, a discussion of indicator value results cannot be conducted with respect to other studies. Our segmentation accuracy results, meanwhile, do correspond with similar studies [100,101], although it is often a challenge to compare segmentation accuracies among GEOBIA case studies due to unique methods, remote sensing imagery, and landscape types.
Although there are no studies to compare this one to directly, we still identified some important limitations. The use of 2D isovist software for visibility calculations means that only quite flat landscapes can be analyzed using the methods shown here. To overcome this, 2D isovist could be replaced by 2D viewshed analysis, which takes the height of a surface model into account; however, as noted by Weitkamp [102], viewshed does not calculate visible space, but visible landscape elements. An algorithm for deriving a viewshed minimal viewpoint set (equivalent to MIS) was proposed by Wang and Dou [103], as well as Shi and Xue [104]. The viewshed approach uses a candidate viewpoint filtering strategy to avoid a rapid increase in the number of viewpoints as the digital surface model resolution increases. Candidate viewpoints with relatively low viewshed contribution rates are iteratively filtered out to derive a minimum set that covers a given viewshed. However, a key limitation of this method is that the VLC of less-open landscapes can be underestimated because the algorithm rejects the candidate viewpoints with a relatively low viewshed contribution rate. Furthermore, depending on the case-study vegetation structure, the 2D viewshed may deviate from reality by up to 45% [105]. This is why, in our case, we preferred to use MIS, despite the fact that it finds a sufficient rather than a minimal viewpoint covering set [56,106]. Even more accurate visibility models can be obtained with the use of 3D-LoS, which uses a similar ray-tracing technique to that used by isovist, or by developing the MIS algorithm to be implemented in 3D isovist software [107,108].
The discussion of the technical limitations of MIS is also relevant for distance decay visibility, view angle, and human eye mimicry. Our methodology uses 360 degrees of observer visibility, assuming that a landscape user (the observer) can look freely in every direction. However, only at the most prominent viewpoints do human observers enjoy 360-degree visibility. The phenomenon of visibility arises in motion, and moving people tend to look predominantly in one direction. Thus, one could argue that an assumption of 360-degree visibility is not particularly natural or realistic for human observers. Agent-based isovist modeling can overcome this limitation. It enables the adjustment of the horizontal view angle (span), as well as the minimum and maximum view range, for simulated observers that are in motion [65]. Therefore, applying an agent-based visibility approach may add additional insights into imageability calculations, expressed not only as a minimum set of vantage points but also as a minimum time needed to overview the landscape. Moreover, because near things are more related than distant things [109], the distance decay effect along an uninterrupted line of sight may also be taken into account during viewpoint creation. To account for this in the field of viewshed modeling, Fisher proposed "fuzzy viewshed" analyses, also called "probable viewshed" analyses [110]. Bartie et al. [111] improved upon this further by taking into account the semi-transparent nature of vegetation; however, neither of these methods were tested in conjunction with MIS analysis.
Another identified limitation with our methodology presented here is imageability indicator value interpretation. Specifically, Vn, VD, and VS did not allow us to capture the process of gradual visual isolation, because these metrics only operate over the entire study area. The increase in Vn and VD values to 2006, followed by their decrease in subsequent years (18% reduction in Vn and VD by 2018) suggests that the MIS algorithm attempts to set up some additional viewpoints to maintain visual connectivity. As it fails (e.g., 2015), two visually separated landscape interiors begin to form VLC branches. However, these results are derived from viewpoint spatial distribution. Similar to the landscape metric biases caused by the land-cover patch scale effect [112], imageability indicators can also be biased by overproduction of secondary viewpoints depending on the adopted scale and LCE smooth factor. In our case, many small, irregularly shaped patches (e.g., shrubs) may have caused the MIS algorithm to overdraw secondary viewpoints.
Our MIS results, supported by the information on land-cover changes from our multi-temporal GEOBIA analysis, show a significant VLC transformation from imageability created by unique landscape openness to imageability as the impression of several landscape rooms enclosed by forest walls. This VLC evolution can be discussed with respect to neuroscience research on the psychological effects of landscape views, providing an example of how imageability indicators can be used beyond merely describing landscape structure. The psychological comfort, pleasure, and stress reduction provided by a long-distance view [58] suggests that similar mental health benefits can be potentially associated with the imageability impression created by primary viewpoints' high field of view. Although our case study site is an uninhabited nature reserve, the calculation of imageability indicators for inhabited areas, as well as associating those indicator values with the quality of a view, may contribute to a better understanding of how to improve subjective well-being [113] concerning urban viewsheds, as well as to better urban planning results. Importantly, further studies on the environmental context in which people are exposed and the mental health implications of those environments [114] could be extended to space-time visual exposure and imageability levels. The landscape users' daily exposure sequences and relative proportions of exposure to primary or secondary viewpoints, recorded by personal tracking devices, could extend knowledge on health-supporting and health-deteriorating factors. Further studies, focusing on people's mobility through a landscape, can address the hypothesis of the mental health benefits of VLC. In order to extend the knowledge of landscape imageability attributes (e.g., landscape layers, color and light, adjacent scenery, visible archetypal elements, landmarks), the artificial expert evaluation tool [115] can be used.
Our findings emphasize the application of GEOBIA to VLC indicator calculations; however, other possible practical applications of the three imageability indicators used here (Vn, VD, VS) can be discussed in the field of landscape quality assessment. The European Landscape Convention [116] obliges signatory countries to assess and monitor landscape quality using landscape indicators derived from quantitative methods [117]. The most common use of landscape quality indicators is based on land-cover patches [118][119][120], which, despite modifiable areal unit biases [21,121], are commonly used in landscape quality studies [122]. Contrary to these studies, we do not evaluate the quality of spatial pattern, but we use land-cover data to conduct visibility analysis (specifically MIS), because it allows us to infer visual landscape structure. Importantly, any quality assessment must be done against pre-established standards and criteria, another area of exploration for landscape imageability analysis. As this pilot study shows, there is strong potential for imageability indicators to find new practical applications, especially in the areas of regional planning, resource management, tourism, or real estate. Another advantage of the analysis shown here is that many aspects of it can be automated. Therefore, it is possible that, in conjunction with Earth observation imagery streams, imageability indicators can be provided as a cloud service similar to other contemporary online land-cover mapping services [123,124].

Conclusions
Despite the technical limitations of using a 2D isovist approach, which can be overcome by 3D isovist software development, we implement in practice, for the first time, the concept of using viewpoints for VLC descriptions. The imageability analysis was itself informed by a multi-temporal GEOBIA land-cover classification process that extends the usefulness of GEOBIA results from ecological indicators and landscape metrics [125] to imageability indicator calculation and derivation. This work contributes to a better understanding of the impact of land-cover change on visual landscape structure dynamics. Our multi-temporal case study reveals considerable changes in case-study VLC resulting from changes in vegetation structure, primarily due to afforestation. Since people identify with landscapes, and since landscapes contribute to psychological well-being [126], VLC is not solely an ecological issue, but also a social one. Thus, lessons learned from how past landscapes changed through time can help to protect and enhance HNV in a more comprehensive way in the future.