Novel Application and Validation of a Method to Assess Visual Impacts of Rotating Wind Turbine Blades Within Woodland Areas

Increasing numbers of wind power plants (WPP) are constructed across the globe to reduce the anthropogenic contribution to global warming. There are, however, concerns on the effects of WPP on human health as well as related effects on wildlife. To address potential effects of WPP in environmental impact assessments, existing models accounting for shadow flickering and noise are widely applied. However, a standardized, yet simple and widely applicable proxy for the visibility of rotating wind turbines in woodland areas was largely lacking up to date. We combined land cover information of forest canopy extracted from orthophotos and airborne laser scanning (LiDAR) data to represent the visibility of rotating wind turbines in five woodland study sites with a high spatial resolution. Performing an in-situ validation in five study areas across Europe which resulted in a unique sample of 1738 independent field observations, we show that our approach adequately predicts from where rotating wind turbine blades are visible within woodlands or not. We thus provide strong evidence, that our approach yields a valuable proxy of the visibility of moving rotor blades with high resolution which in turn can be applied in environmental impact assessments of WPP within woodlands worldwide.


Introduction
Addressing global warming and related counteractions, a number of international political agreements aimed to promote the use and production of renewable energy in form of wind energy, photovoltaics or other techniques within the last decades (UNFCCC 2015). As a consequence, wind power plants (WPP) have been increasingly constructed across the globe (Renewable Energy Network 2018). There are, however, increasing concerns regarding the effects of WPP on human welfare (Merlin et al. 2015) and on free living animals (Kuvlesky et al. 2007).
Assessing and analyzing environmental impacts of planned WPP on humans typically account for acoustic effects (noise emission), visual effects (e.g. shadow flickering, reflecting or blinking) or risks in form of ice throw on humans or human infrastructures (e.g. Kuvlesky et al. 2007;Schäffer et al. 2019). Residential health aspects or psychosocial stress as well as landscape aesthetics are further main issues that have been addressed in past studies and assessments (Bakker et al. 2012;de Vries et al. 2012). Such impact assessments of WPP typically focus on human wellbeing, including the number of turbines, distances to residential or recreational functions, observers' height, weather/atmospheric conditions or forest cover acting as screen (Manchado et al. 2019;Möller 2006;Rodrigues et al. 2010;Torres et al. 2009;Tsoutsos et al. 2009; for an overview see also Jombach et al. 2010). However, comprehensive environmental impact assessments should also address effects of WPP within woodlands on wildlife.
Effects of WPP on wildlife include the collision of bats and birds with moving wind turbine blades, barotraumas (De Lucas and Perrow 2017;Drewitt and Langston 2008) as well as behavioural adaptations due to WPP presence (Coppes et al. 2020;Hötker 2017;Rabin et al. 2006). These adaptations include changes in birds' vocalization (Whalen et al. 2019;Zwart et al. 2016) and the avoidance of areas close to the turbines (Coppes et al. 2020;Hötker 2017). It is assumed that noise generated by WPP causes changes in displaying (singing) behaviour of male birds (Szymański et al. 2017), but animals seem also to react to visual cues (D'Angelo et al. 2008), e.g. shadow flickering or moving blades of WPP.
While several studies focus on effects of WPP on wildlife behaviour in flat, open landscapes with low amounts of vegetation (Hoover and Morrison 2005;Millon et al. 2015;Smith et al. 2020), large knowledge gaps remain for forested areas, although there is an expanding pressure on these habitats (Bunzel et al. 2019). However, some rare studies provided evidence, that the movement of rotating turbine blades yielded behavioural effects on forest-dwelling wildlife species (Coppes et al. 2019;Skarin et al. 2018). Movement of rotor blades in the view of animals' vision might trigger an anti-predator response, as prey animals react to movement in their sight (Heesy 2004). As a consequence, animals avoid such objects to reduce the risk of predation. Thus, home ranges of freely ranging female reindeers shifted after WPP construction and switched to areas where wind turbines were "out of sight" (Skarin et al. 2018). Thereby, "out of sight" was defined based on topography at a target height of 110 cm above ground, including information on vegetation coverage. Such studies are generally known as visual impact assessments.
Noise and shadow flickering of WPP can be modelled for WPP using specialized software (e.g. windPRO-EMD International). One widely applied approach for visual impact assessments of WPP are viewshed analyses (cf. Lagner et al. 2018). This method draws lines of sight between WPP and selected locations in the landscape (Möller 2006). As large obstacles affect visibility of the WPP (Burrough et al. 2015), this analysis generally includes the topography, buildings and in some cases the average height of forests. Especially in forested landscapes, the structure of a forest (i.e. open vs. dense canopy cover and tree height) is likely to affect from where a WPP is visible or not. Approaches accounting for such woodland characteristics were only scarcely published yet (Skarin et al. 2018) or are entirely missing in such visual impact assessments. However, to obtain reliable estimates of visual cues of WPP on forest-or ground-dwelling animals, it is crucial to determine an accurate assessment of the visibility of rotating turbine blades within woodlands.
Facing this lack, we developed a high-resolution visibility assessment of rotating turbine blades within forest-covered landscapes. We combined land cover information on forest canopy extracted from orthophotos and airborne laser scanning (LiDAR) data to assess the visibility of rotating turbine blades within forested landscapes on five study sites with high spatial resolution. Thereby, we aimed at the construction of a valid visibility assessment, which can be derived from earth observation data rather than field records and which can be taken as proxy for the visibility of rotating turbine blades in forest landscapes. Our approach allows for broad application by practitioners to ensure its consideration in environmental impact assessments in forest-covered landscapes. Conformingly, we used broadly available data without the need for extensive preparative procedures. This ensures practicability of the assessment and thus enables subsequent consideration of the resulting proxy in environmental impact assessments. As such visual impact assessments of moving rotor blades within forested areas have not been considered in environmental impact assessment up to now, our approach closes this lack and provides a further tool for decision making. In order to evaluate our visibility assessment for its usability, we developed a comprehensive validation approach. We used a unique, large dataset of independent field records of the visibility of rotating turbine blades, including different forest sites across Europe. This validation was done by means of Generalized Linear Models (GLM), where we contrasted visibility-predicted by our assessment approach-with levels of visibility recorded in nature. To our knowledge, both aspects, i.e. the highresolution assessment of the visibility of rotating turbine blades within woodlands and the comprehensive statistical validation of the assessment based on cross-European field

Study Sites
We used five cross-European study areas with active WPP comprising study sites in the Eastern Alps (Austria), on central European low altitude mountain ranges (Germany) and Central Sweden (see Table 1, Fig. 1

Visibility Assessment (VA)
Orthophotos for the land cover classification and topographic data based on LiDAR technology were supplied by the Styrian provincial government ("Land Steiermark, A17 Statistik und Geoinformation, GIS", Austrian study area), the regional authority for geoinformation and land development Baden-Württemberg ("Landesamt für Geoinformation und Landentwicklung Baden-Württemberg", German study areas) and the Swedish mapping, cadastral and land registration authority ("Lantmäteriet", Swedish study area) ( Table 2). We used the digital surface models (DSM) and digital terrain models (DTM), prepared and provided by the Fig. 2 Workflow of data processing, field sampling and statistical validation within our study provincial governments or regional authorities, respectively, as this would allow for a congruent processing in practical application contexts in future. For the Austrian Province of Styria and for Sweden, topographic data are provided as open data.
To address effects of forest cover and related vegetation structures on the visibility of rotating turbine blades, we first defined a processing radius of 3 km around each turbine and then derived the height of the vegetation for all study areas. We limited the processing radius to 3 km to reduce calculation efforts and to focus on an area of presumably higher biological relevance (i.e. decreasing perceptibility of rotating turbine blades with increasing distance). Using the highresolution topographic data based on LiDAR technology,   Fig. 2). Additionally, we classified false-colour orthophotos in a pixelwise manner by thresholding the bands. This straightforward classification approach was sufficient as the dominant land cover was forest and the purpose of the classification was to eliminate forest from non-forest land cover. The thresholding had to be done for each study site separately.
In a second step, we only considered vegetation higher than 2 m in our analyses to further eliminate the phenological highly variable ground vegetation. Therefore, we could separate this stratum, which does not significantly impact the visibility of rotating turbine blades for medium-to largesized forest-dwelling wildlife species, from forest cover in a narrow sense. We then added only the resulting forest vegetation height to the DTMs (i.e. elevation and topography of the study areas) to create a raster dataset (1 × 1 m resolution).
Finally, we applied a viewshed analysis per turbine using the tool "Viewshed" in ArcGIS 10.5.2 with the following input parameters: the geographical location of each single turbine, the total height of the turbine and our raster including data on elevation, topography and forest vegetation. To create and validate a standardized, yet simple and widely applicable proxy for the visibility of rotating wind turbines in woodland areas, we calculated visibility at ground level per turbine. The output of the viewshed data indicates per raster cell (1 × 1 m resolution), whether the upper tip of the rotating turbine blades is visible or not (binary output) from this cell (see Fig. 3). Thus, instead of modeling the visibility of the entire rotor swept area (e.g. by combining several viewshed analyses), we created a model which is readily applicable and easy to interpret (Fig. 4).

Field-Based Visibility Records (VR)
For a validation of our VA, field work was performed between 1 July 2016 and 31 October 2016. We recorded visibility of rotating turbine blades within woodlands via field To verify, whether our VA of rotating turbine blades would hold as authoritative proxies for impact assessment for medium-to large-sized forest-dwelling animals, we recorded the visibility of moving blades at two different observer heights above ground (explained below). The fieldbased visibility records were based on several, subsequent steps: first, we compiled a schematic representation of three visibility classes (not visible/partially visible/fully visible), which are highlighted and described in Table 3. Prior to the field records, we intensively trained all field observers to calibrate field records and to ensure consistent assessments throughout the study. During the field records, observers compared and recalibrated their assessment in several day intervals by comparing and adjusting their personal assignment in accordance with other observers. This process ensured a consistent recording throughout the study.
To account for positioning inaccuracy of GPS locations, we based our visibility assessment on a 5-m contour line transect-average per sampling plot. Field records were only done on days with optimal weather conditions (no fog, no rain) representing time of maximum visibility and when rotor blades were rotating/active. We used an observer height of (1) 50 cm above ground, which simulates the view angle of different medium-sized ground dwelling animal species in Central European forests (e.g. red fox Vulpes, capercaillie Tetrao urogallus) and (2) 150 cm above ground, which simulates the view angle of large ground-dwelling animals (e.g. red deer Cervus elaphus) and also represents human perception (cf. Lagner et al. 2018). The validation of the visibility assessments at these two observer heights should indicate, whether the general viewshed analyses at ground level would provide an authoritative proxy of the visibility for a wide range of ground-dwelling wildlife species of differing body size. To account for a potential impact of different environmental drivers on the visibility of rotating turbine blades within woody landscapes and to address potential drivers of misclassification, we also recorded a set of explanatory variables during the field records (Table 4).

Contrasting Assessed Visibility vs. Field Records-Statistical Validation
We validated our visibility assessment (VA) of rotating turbine blades by contrasting the modelled visibility against the results of our independent field-based visibility records (VR). VA was determined at different spatial scales (raster cell resolution 1 × 1 m, within a 10 and 20 m radius), whereas VR was determined at two different observer heights (50 cm and 150 cm above ground). Field validation was only done within woodland areas. Based on the binary VA output per grid cell (not visible/visible), we calculated the spatial percentage of cells, where rotating turbine blades were predicted to be visible. This percentage of cells with visibility was derived for a 10 m and 20 m radius around each field sampling point, using a circular moving window. This further improved our model validation, as the grid points of the field records were approached using handheld GPS, which are known to hold potential GPS-inaccuracies of up to several meters, while the grid cell size for the viewshed analyses (not visible/visible) was 1 × 1 m. Considering the ordinal and binary scale of our VR output (i.e. three classes vs. two classes of visibility, see Table 5) and the different scale levels of potential drivers of visibility in the field (see Table 5), we chose GLMs (Generalized Linear Models) and GLMMs (Generalized Linear Mixed Models) for the statistical validation. Logistic regression models describe the relationship between a categorical response variable and at least one explanatory variable (Anderson 1984;Hosmer and Lemeshow 2000). In case that the response variable is dichotomous, a binary logistic regression ca be run. A special case of multinomial logistic regression models, where the response variable shows more than two categories (being polytomous), is the ordinal logistic regression. For this type of logistic regression, the response variable shows at least three ordered levels. For proportional odds ordinal logistic regression models, which we used in our study, the proportional odds Occurring/not occurring (0/1) at sample plot level Canopy cover % of the ground covered by the crown of trees assumption has to be satisfied. According to this assumption, the slope between each pair of levels of the response variable is assumed to be the same for each independent variable. Thus, the model comprises only one set of coefficients for the independent variables, but different intercepts per level of the response variable, respectively. In general, logistic regression models hold the advantage that they are robust both in terms of underlying distributions and in terms of scale levels of input data (Agresti 1989). The statistical validation of the visibility assessment with GLMs was based on two steps (see also Table 5): 1. We first calculated proportional odds ordinal logistic regression models (OLR) per study site for each observer height (50 cm and 150 cm above ground) separately, using the three visibility classes (not visible/partially visible/fully visible) of the field-based visibility records as response variable and the calculated percentage of visibility as explanatory variable (VA, within a 10 m or within a 20 m radius, respectively). For the study site Alexanderschanze, we did not calculate an OLR, as the sample size per visibility class was too small to allow for site-specific OLR runs (Harrell 2001). We further fitted models with additional environmental variables (e.g. wind conditions or backlighting, for details see Tables 4 and 7) to test the explanatory value for the visibility records (VR). Prior to these analyses, we plotted the explanatory variables against the visibility assessment (VA) classes and run Brants tests (Brant 1990) to ensure, that the explanatory variables would satisfy the proportional odds assumption, being one precondition for running OLRs (McCullagh 1980;Harrell 2015). We also compared the effect coefficients of environmental explanatory variables based on binary logistic regressions of two related levels of the factorial response variable, respectively (i.e. not visible vs. partially visible and partially visible vs. fully visible) to control for the proportional odds assumption. As the environmental explanatory variable "canopy cover of trees" did not meet the requirement of proportional odds, it was discarded from OLR calculations. 2. Accounting for the fuzziness of visibility records (VR) in the field (i.e. the field workers' assignment of the visibility of rotor blades to one of three distinct classes of visibility-not visible/partially visible/fully visible), we also ran mixed effects logistic regression models with a binary outcome variable (binary logistic regression models, BLR) considering only two classes of visibility (not visible vs. visible) at a given study plot. Thereby, we used the binary visibility records (VR) (at 50 and 150 cm above ground, respectively) as response variable, the visibility assessment (VA, within a 10 m or 20 m radius) as well as additional environmental varia- Table 5 Description of the different validation models (model type, related input data and geographic resolution are indicated with • or ○,

respectively)
Environmental variables, which were tested, but did not enter the final models, are indicated as ○.
Environmental variables, which were not tested within the ORLs (not satisfying the proportional odds assumption) are indicated as ---Vali-dation models Visibility-field records Visibility assessment-output bles (e.g. wind conditions, backlighting or canopy cover of trees) as explanatory variables. Accounting for the differing sample sizes per study sites (Table 1), we considered the study site as random intercept (after checking the random intercept assumption via boxplots per study site).
To analyze potential causes for misclassifications of the VA (prediction of the VA does not correspond to the VR), we further ran BLRs. We used correctness of classification (i.e. correct assignment of plots with or without visibility of rotor blades) of the VA as binary response variable and different environmental variables (e.g. wind conditions, backlighting and canopy cover of trees) as explanatory variables. Again, we considered the study site as random intercept, accounting for differing sample sizes per study site.
The models were calculated using the packages MASS (Venables and Ripley 2002) for the OLR and lme4 (Bates et al. 2015) for the BLR within the R 3.5.2 (R core team) environment. The Brants test was calculated using the brant package (Brant 1990).

Results
In total, we validated the output of our visibility assessment with 1738 independent direct personal recordings of visibility in the field. Both, across all study sites and for each single study site, the percentage of visible grid cells (according to the model) within a 10 m and a 20 m radius increased with increasing visibility levels recorded in the field (not visible/partially visible/fully visible) (for an overview across all study sites see Fig. 5). This applies both for observer heights of 50 cm and 150 cm above ground (Fig. 5). At plot center level (1 m 2 ) and for both observer heights above ground (50 cm and 150 cm), the rate of correct classification (= both status "visible" and "not visible" according to the model and to the field records) was 69% across all study sites (n total_plots = 1738). Correct positive classification rates (visibility = correctly classified) was 67 and 69% and correct negative classification rates (non-visibility = correctly classified) was 70 and 68% for 50 cm and 150 cm observer heights, respectively.
At study site level, OLRs consistently yielded a significant explanatory effect of the visibility assessment (VA) on the visibility classes (not visible/partially visible/fully visible), recorded in the field (VR). This applies both for the 10 m and the 20 m radius and for observer heights of 50 cm and 150 cm above ground (see Table 6, "Appendix"). Apart from the VA, only the wind condition appeared as significant environmental explanatory variable in rare cases. Neither sun conditions (backlighting), nor the occurrence of forest edges significantly contributed as explanatory variables to the OLRs per study site. We, thus, reduced the OLRs per study site to the intercepts and the explanatory variable VA.
Across all study sites, the BLRs showed significant effects of the VA (10 m and 20 m radius) on the visibility classes of VR, together with effects of the environmental explanatory variables canopy cover of trees and wind conditions (see Table 7, "Appendix"). The canopy cover of trees negatively affected visibility of rotor blades in the field, whereas increasing wind stringently showed a positive effect on VR. Sun conditions (backlighting) and the occurrence of forest edges did not significantly contribute to the BLRs across all study sites.
The BLRs of classification correctness also showed significant effects of the environmental explanatory variables canopy cover of trees and wind conditions on the classification correctness (i.e. correct classification of both visibility and non-visibility; see Table 8, "Appendix") at both observer heights. Thereby, both increasing canopy cover of trees and increasing wind negatively affected accordance between VA and VR. Again, the explanatory variables of sun conditions (backlighting) and the occurrence of forest edges did not enter the final BLRs.

Discussion
The validation steps of this study showed, that our visibility assessment provides an authoritative proxy for the visibility of rotating turbine blades within landscapes that are covered by forests. This applies both for different levels of visibility of rotating turbine blades (not visible/fully visible/partially visible) and for a binary assessment (visible vs. not visible). Thus, we offer a valid proxy that had neither been derived nor statistically validated for forested landscape patches until now.
Independent field observations of the visibility and related model validations for woodland areas are scarcely available in literature (cf. Wróżyński et al. 2016). In all cases within our study, the visibility assessment proved to be a significant independent variable within the different validation models. Consequently, high resolution land surface data generated via LiDAR technology, can adequately depict forest and vegetation structures which lower the visibility of rotating wind turbine blades apart from pure terrain features. Using our experimental approach, we proved that natural visual obstructions (e.g. forest vegetation cover) can be considered in spatial visibility analyses. Thus, neither time consuming, expensive field records, nor direct addressing of residents with questionnaires (Pedersen et al. 2009;Tsoutsos et al. 2009) or related methods (de Vries et al. 2012) have to be conducted.
The output of our visibility assessment can be used to implement the visibility of rotating turbine blades adequately in environmental impact assessments in woody landscapes, which was lacking so far (e.g. Rafiee et al. 2018;Rodrigues et al. 2010) and to base conclusions on potential effects on the visibility assessment output. This might address both, human welfare perspectives as well as potential effects of WPP on forest-dwelling wildlife (see Coppes et al. 2019), constituting a potentially important factor in regional assessments of visual impacts and in related planning processes (cf. Möller 2006). As distinctly different "observer heights" might be relevant for considering effects of rotating turbine blades on different medium-to large-sized ground-dwelling wildlife-species, the validation of the ground-level VA at different observer heights was an important issue in this study. Thereby, we could prove, that even ground-level VA provides an authoritative, yet simple and widely applicable proxy for the visibility of rotating wind turbines in woodland areas.
However, despite of a high rate of correct classification (and thus a successful model validation), some uncertainties in the purely remote-based assessment of visibility are remaining. All final binary logistic regression models (at different observer heights and for different plot radii), either proving relations of the visibility assessment output to the visibility records in the field or exploring driving factors of misclassification, comprised the variables canopy cover of trees and wind condition. For the first suite of models, relating the output of the visibility assessment to the field records, increasing canopy cover within forests negatively impacted visibility in the field, whereas increasing wind had the opposite effect. Obviously, the assignment of sample plots to the category "visible" or "not visible" becomes fuzzier in case of higher canopy cover. This applies for both validation heights and for the two plot radii. In contrast, fast blade motion during strong wind conditions might attract views. The referring effect coefficients in the validation models all show a very similar order of magnitude.
On the other hand, both canopy cover and wind affected the classification correctness, with a lower probability of correct classification in case of increasing forest cover and increasing wind. While visibility on sparsely covered plots might be easily recorded, increasing forest cover could have caused difficulties for the field workers to determine a distinct visibility class. Another possible cause for the observed misclassifications of the visibility could be attributed to the fact, that our visibility assessment only addressed the visibility of the tip of the WPP and not of the entire rotor swept area. Therefore, the upper tip of the turbine blades might not be visible as per assessment, but a large part of the rotor swept might be visible in the field recording, for example at forest edges. However, forest edges did not remain as significant driver in any final validation model.
We finally conclude, that our assessment of the visibility of rotating turbine blades in woodlands can adequately estimate from where moving blades are visible or not. Based on a unique field validation of the assessment in different topographical conditions throughout Europe, we assume that our approach can be applied under different forest conditions. We thus provide an authoritative additional approach to represent potential impacts of WPP in woodland areas on both humans and wildlife in terms of visibility of rotating turbine blades. This information can consequently be used for evaluations of existing as well as planned WPP in woodland areas.

Appendix
See Tables 6, 7 and 8. Table 6 Final ordinal logistic regression models per study site at the two different observer heights (50 cm and 150 cm above ground); observer heights are indicated with "_50cm" and "_150cm" in the model acronyms; study sites are indicated in the model acronyms with Hochpürschtling, HP; Pretul, PR; Hornisgrinde, HG; Sweden, SWE; categories of the response variable (visibility of rotor blades during field records) are: 1, not visible; 2, partially visible; 3, fully visible; the output of the visibility model VM is indicated as "VM_10m" and "VM_20m", indicating the percentage of 1 m 2 grid cells within a 10 or 20 radius, where the visibility model yielded visibility of the rotor blades; intercepts are depicted per category of the response variable, relating the levels 1 (partial visibility) and 2 (full visibility) to the reference category 0 (no visibility) Model

Compliance with ethical standards
Conflict of interest All authors declare that they have no financial or personal relationships with other people or organizations that could inappropriately influence their work and declare that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, Exp(B) denotes the effect coefficient (= odds ratio), SE the standard error associated with the effect coefficient; t value is the parameter estimate (aka coefficient) divided by its standard error and p is the p-value used in testing the null hypothesis that the coefficient (parameter) is 0  Table 7 Final binary logistic regression (BLR) models across all study sites at the two different observer heights (50 cm and 150 cm above ground); observer heights are indicated with "_50cm" and "_150cm" in the model acronyms; study sites considered as random intercepts; categories of the response variable (visibility of rotor blades during field records) are: 1, not visible; 2, visible; the output of the visibility model VM is indicated as "VM_10m" and "VM_20m", indicating the percentage of 1 m 2 grid cells within a 10 or 20 radius, where the visibility model yielded visibility of the rotor blades Wind conditions are indicated as "wind" with three distinct wind classes, i.e. wind 1, Bft 0-1; wind 2, Bft 2-3 and wind 3, Bft 4-5. Exp(B) denotes the effect coefficient (= odds ratio), SE the standard error associated with the effect coefficient; t value is the parameter estimate (aka coefficient) divided by its standard error and p is the p-value used in testing the null hypothesis that the coefficient (parameter) is 0  Table 8 Final binary logistic regression (BLR) models of classification correctness (across all study sites), calculated per sample plot center at the two different observer heights (50 cm and 150 cm above ground); observer heights are indicated with "_50cm" and "_150cm" in the model acronyms; study sites considered as random intercepts; categories of the response variable (correctness of classification) are: 0 = not correctly classified (i.e. either visibility or non-visibility), 1 = correctly classified (i.e. both visibility and non-visibility). Wind conditions are indicated as "wind" with three distinct wind classes, i.e. wind 1, Bft 0-1; wind 2, Bft 2-3 and wind 3, Bft 4-5 Exp(B) denotes the effect coefficient (= odds ratio), SE the standard error associated with the effect coefficient; t value is the parameter estimate (aka coefficient) divided by its standard error and p is the p-value used in testing the null hypothesis that the coefficient (parameter) is 0 adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.