Spatial Analysis of Ecological Risk in a Coastal Dune Landscape Using High Resolution Aerial Photography

This study was conducted in Parangtritis coastal dune, a landscape located in vicinity to Parangtritis village in the southern part of the Special Region of Yogyakarta of Java Island, Indonesia (Fig. 1). As for the absolute location, it spans from 110°17’0” E, 8°0’30” S to 110°20’30” E, 8°1’30” S. The dune area is bordering with a wide beach and Indian Ocean from the south, the Opak River from the west, a highland with a high cliff from the east, and a partly forested and settled agricultural land from the north.

Fig 1.

The Parangtritis coastal dune area is 412.8 ha, consisting of the supporting zone in the east (176.4 ha), the centre or core zone (141.1 ha) and the restricted zone in the west (95.3 ha). It is a unique landscape as it is the first-discovered and one of the last not urbanised coastal dunes containing a barchan morphology in Southeast Asia (Sunarto 2014, Nehren et al. 2016). This barchan has been formed in the tropical-humid climate zone, specifically in the coastal area, instead of in the arid zone where it is usually formed (Sunarto et al. 2010).

Its morphology consists of small coastal dunes with wind migrating towards southeast barchan forms, and then entering the coastal agricultural and settled (inhabited) land. The beach has tide and wave which influenced/inducted the formation of the coastal dunes. The tide and wave move and reach the beach covered by sand which has 40–80 m width. The barchan forms are developed due to wind-driven transport from the beach to the land with southeast direction. The height of the barchan nearly reached 11 m in 2011 but decreased 2 m in 2020.

This ecosystem is an almost natural habitat characterised with a moderate climate (D) based on Schmidt–Ferguson classification of the meteorological statistics (Putri 2008), an annual rainfall of 1000–2000 mm (Services of Public Works, Housing, Energy and Mineral Resources for the Special Region of Yogyakarta 2014), an annual temperature of 26–27°C (Malawani 2014) and an average wind speed of 5.3–9.2 m s⁻¹ (Services of Public Works, Housing, Energy and Mineral Resources for the Special Region of Yogyakarta 2014). The word almost presents before the natural habitat because other than the native flora and fauna, there are also introduced flora and fauna on the sand dunes. Moreover, human intervention could be seen by their activity and the distribution of non-bare land cover observed rapidly from the aerial photography, mainly within the supporting zone and restricted zone.

The study area is a habitat for some native desert-like flora, but coastal and tropical climate zones are characterised by torn-like small-leaf and long-root flora. Other than native flora, there are also introduced flora such as Eupatorium inulifolium, Acacia mangium, Swietenia mahagoni, Anacardium occidentale, Gliricidia sepium and Casuarina equisetifolia. The non-native flora is confirmed to be a disturbance agent for the coastal dune because it has a wetter condition that will affect the creation of a new ecosystem that is not suitable with the original coastal dune condition (Services of Environment for the Special Region of Yogyakarta 2017). Besides, the introduced flora tend to indirectly change the dune morphology through the change in process intensity.

The Parangtritis coastal dune materials come from Merapi volcano, transported through rivers to meet the Opak River as the mouth before entering the Indian Ocean. The wave and current move the sand from the ocean to be deposited on the beach. After that, the wind transports the dry sand from the beach to the land. The wind that brings the dry sand will slowly deposit the sand on the land, forming dunes, when the energy becomes lower and the movement is trapped by the Batur–Agung Hills (Fig. 2). Since 1985, Fakhruddin et al. (2010) reported that the Kretek Bridge over the Opak River has stimulated the development of infrastructure in Parangtritis to support its function as a tourism location, including the settlement. The material supply from the Opak River reduces with time. It is also aggravated by the buildings and Casuarina forest along the beach, which become an obstacle for the wind to bring the sand to the land.

Fig 2.

A landscape ecological units (LEU) is a unit in analysing an ecological system or geosystem, which consists of geology, relief and natural conditions creating the materials within the land (Burel, Baudry 2003). The LEU can be determined based on the terrain landscape, non-terrain landscape (such as LULC) or a combination of the two (Chattaraj et al. 2018). This study assumed Parangtritis coastal dune as one terrain landscape coastal dune, while the next level landscape is characterised by the LULC for detailing the LEU.

The landscape ecology has a relation with anthropogenic geomorphology (Csorba 2010). Therefore, the anthropogenic impact to the landscape containing human intervention proposed by Szabó (2010) can be considered a determination method for the LEU. In this study, LEUs made by factoring in human intervention are called landscape class, while the one designed based on LULC is termed landscape sub-class. Small format aerial photographs (SFAPs) in 2011 and 2020 with spatial resolutions of 30 cm and 10 cm, respectively, were visually interpreted to identify the sub-class. The minimum mapping unit applied to the image interpretation for a 1:2000 map is about 158.76 m². A map of the coastal dune obtained from the Parangtritis Geomaritime Science Park was also consulted in the SFAP interpretation, with some adjustment to the national classification system. Each identified LULC was complemented with human intervention types proposed by Szabó (2010). The interpretation results, i.e., LEU maps, were later verified through a fieldwork and then re-interpreted and revised.

A landscape ecological risk was assessed by calculating the landscape disturbance and landscape fragility (Xie et al. 2013) using the landscape class and sub-class scenarios with the pivot tables in a spreadsheet program that were integrated with a set of tabular geographic information system. Landscape disturbance (S_i) indicates the degree of interference in different landscapes and for this reason considers existing human activities. It was calculated based on the detected landscape pattern, including fragmentation, isolation or dominance. The landscape fragmentation index (C_i) measures the complexity of a segmented landscape’s spatial structure and the degree of disruption due to human activities. The landscape isolation index (N_i) reflects the degree of a patch distribution for a certain landscape type. The landscape dominance index (D_i) explains the degree of the patch’s influence on the landscape pattern formation and changes. The equations used to calculate these indices are presented in Table 1.

Landscape fragility (F_i) in this study was assessed using an analytical hierarchy process (AHP) (Liu et al. 2020), which was performed using the spatial multicriteria tool in ILWIS 3.4, i.e., pairwise comparison. The landscape fragility index is the priority value belonging to each class, which is accepted because the consistency ratio is 0.0754 (<0.1). The fragility index of each class was applied to the sub-class (Xie et al. 2013, Shi et al. 2015, Liu et al. 2020) because the AHP runs properly only with a limited number of pair comparisons (Ozdemir 2005). Further explanation is available in a previous study (Listyaningrum et al. in press), with the summary presented in Table 2.

The ecological risk index (ERI) of the identified landscape types was calculated using the equation in Table 3, as formulated by Zhang et al. (2013). The assessment unit of the ERI was a grid of equal size whose value was represented by a centroid (Fig. 3). The grid is suitable for assessment in a small area, such as the coastal zone observed in this study (Liu et al. 2020). The grid was 63 × 63 m in size, which is the minimum area mappable or observable in a 1:10,000 assessment.

Fig 3.

Spatial interpolation enables ecological modelling, which transforms vector data, including centroid, into raster data with continuous values. Some aspects to consider while selecting an interpolation technique are the sample data, type of surfaces to be generated and tolerance of the estimated errors (Kumar, Sinha 2018). This study selected Kriging (indicator type) because it has the least root-mean-square error (RMSE) among interpolation techniques (Requia et al. 2019), which satisfies the spatial distribution of the ERI value in the study area.

This study analysed the LER using global and local spatial autocorrelations (Xie et al. 2013, Zhang et al. 2018). The former measures the degree of global correlation and disparity between ecological phenomena using the common indicator Moran’s I, while the latter uses a local indicator of spatial association (LISA) to illustrate the spatial distribution of local heterogeneity and determine the degree of the spatial disparity between a region and its surroundings. The LISA detects any significant spatial clustering created by a region and its peripheries (Anselin 1995). The formulas used in the global and local spatial autocorrelations are shown in Table 4.

The spatial distribution of LULC in the study area is illustrated in Fig. 4. There were 34 and 37 sub-classes of LEUs in 2011 and 2020, respectively, corresponding to the detailed classification in Liu et al. (2020). The study area was mainly covered by bare land and sparse shrubs. In 2011, the bare land was more expansive (29.84% of the total area) than the sparse shrubs (19.89%), making it the landscape matrix for this year. However, in 2020, it had decreased significantly and became the second-largest LULC (19.78%) after the sparse shrubs (25.68%).

Fig 4.

Besides the area metrics, aggregation metrics such as the number of patches were also considered in the LER calculation. The number of patches in the study area varied from 1 to 278 in 2011 and from 1 to 338 in 2020. In 2011, the bare land had the highest number of patches (278), followed by sparse shrubs with 225 patches, whereas in 2020, sparse shrubs had the most patches (338), followed by bare land, which amounted to 253. Figure 4A, B also show that the landscape became more diverse significantly including (1) the addition of new LULCs such as fisheries and main road to the supporting zone; (2) in the core zone, the bare land decreased in size but the forest area increased in size; and (3) the bare land in the restricted zone also decreased, and the one in the southeast was transformed into fisheries and tourist sites.

The spatial distribution of human interventions is presented in Fig. 5. Parangtritis coastal dune was grouped into natural landscape and anthropogenic landscape; the latter was subdivided into eight classes. In the course of 10 years from 2011 to 2020, natural landscape dominated the dune, covering up to 64.10% and 51.25% of the entire area, respectively.

Fig 5.

The number of patches in the study area was in the range of 1–118 in 2011 and 3–246 in 2020. Natural landscape comprised the highest number of patches in 2011, which later lowered to the second-highest in 2020 after agrogenic landscape. Figure 5 illustrates that tree agrogenic was the main landscape class in the supporting zone, natural in the core zone and urbanogenic and tourism-sports in the restricted zone. In 2010, industrogenic was found in a small proportion in the supporting zone before expanding and sprawling to the restricted zone in 2020. Similarly, the landscape formed by info-telecommunication interventions showed a similar pattern of change to industrogenic. Anthropogenic interventions such as narrow paths called corridors formed the traffic landscape, which was distributed everywhere. Water management was a minor component of the core zone and restricted zone, which, like the traffic, formed corridors too.

The LER assessment produced ERI values ranging from 0.0 to 1.0, which are spatially distributed in Figure 6. To understand the ecological risk level, these values were equally divided into very low (ERI ≤ 0.2), low (0.2 ≤ ERI ≤ 0.4), moderate (0.4 ≤ ERI ≤ 0.6), high (0.6 ≤ ERI ≤ 0.8) and very high (0.8 ≤ ERI ≤ 1.0) (Gong et al. 2015). Figure 7 shows that all these LER levels changed insignificantly in their spatial distribution, except for the very low risk in the 2011 landscape class (LEU) scenario.

Fig 6.

Fig 7.

Table 5 shows that for each scenario of landscape classes and sub-classes, Moran’s I was higher than zero, ranging from 0.4112 to 0.6342. It indicated that within the 63 × 63-m grid, there was a positive spatial autocorrelation between the LER indices. A relatively high Moran’s I also suggested clustering of LER in the study area, except for the 2011 sub-class scenario whose Moran’s I was 0.4112 (< 0.5). This may have caused the spatial distribution of the ERI values in the 2011 subclass to be different from those in other scenarios, specifically in the core zone (Fig. 6A). Overall, both class and sub-class scenarios showed that Moran’s I increased in the years 2011–2020, signifying an increasingly stronger global spatial autocorrelation. The minimum Z_score in the study area was 18.0872 (> 1.96). As a result, at a significance level of 0.05, the null hypothesis was rejected. In other words, the values of the spatial attributes have spatial autocorrelation.

The local spatial autocorrelation of the LER in the study area is presented in four LISA cluster maps in Figure 8, which enable the analyses of Figure 6 in a cluster view and the significance of a unit’s LER level on its neighbours. The autocorrelation process yielded high–high and low–low clusters, high–low and low–high outliers and not significant in 63 × 63 m grids. Each scenario has specific LEUs and patterns depending on the cluster formed.

Fig 8.

Image interpretation elements in remote sensing, landscape structure in the landscape ecology, the LULC classification system and the human intervention classification system have been considered in defining the LULC types in this study. Referring to Lillesand et al. (2015), the eight keys of image interpretation are size, shape, shadow, tone, colour, texture, pattern and association. Because the landscape of the Parangtritis coastal dune is characterised by linear objects such as corridors and non-linear objects such as matrixes and patches (Forman 1995, Burel, Baudry 2003), shape is the first element considered in interpretation, instead of colour/tone. The LULC classification system has been developed using the national standard (National Standardisation Agency 2010, Geospatial Information Agency 2016), with some modifications made to acquire spatial data and attributes of sub-classes in 2011 and 2020.

In terms of LULC classification details, monitoring environmental elements with SFAP using the Szabó (2010) classification system simplifies the abiotic, biotic and cultural elements of LULC into two: abiotic and cultural. While non-built-up areas belong to the abiotic element, built-up areas are an example of the cultural element. Because fauna cannot be expressed in SFAP, the biotic element only comprises flora or vegetation. The flora in the study area are divided into native vegetation (abiotic) and introduced vegetation (cultural).

Szabó (2010) classified human interventions into montanogenic, industrogenic, urbanogenic, traffic, water management, agrogenic, warfare, tourism and sports (in this study, tourism and sports are combined into tourism-sports). Montanogenic and warfare interventions are omitted because they do not exist in the study area. Overall, only human interventions that cause morphological changes like excavation, plantation or accumulation are observed. The derived LULC map has been used to investigate these alterations (Adzima et al. 2020).

The research factor is natural landscape because not the entire coastal dune is influenced by human interventions. The natural landscape consists of bare/ vacant land, grassland, sparse shrubs, footpaths, litter accumulation and natural river. Bare land is the priority ecological unit in the landscape because it provides natural habitats. Grassland and sparse shrubs are grown with native vegetation species that do not disturb the coastal dune formation (Services of Environment for the Special Region of Yogyakarta 2017). Footpaths and litter accumulation also do not significantly influence the morphology and the natural process occurring in the coastal dune. Natural river is distinguished from seasonal river because the former occurs naturally, whereas the latter is a construction installed on river banks with water management functions.

The class agrogenic has been modified and subdivided into agrogenic (general agrogenic) and tree agrogenic to accommodate varying species of introduced vegetation in the study area. Tree agrogenic is separated from agrogenic because the local multiple-species and Casuarina forests have a relatively high surface (canopies), which together with green paths and dense shrubs create windbreakers on the coastline (Syahbudin et al. 2013a, b) and potentially block sand-carrying wind as the primary agent in coastal dune formation. Moreover, some vegetation species are not native to the Parangtritis coastal dune ecosystem such as A. mangium, S.mahagoni and A. occidentale in the multiple-species forest; C. equisetifolia in Casuarina forests and green paths; and G. sepium in the green paths and dense shrubs (Oktavianto, Handayani 2017, Services of Environment for the Special Region of Yogyakarta 2017, Widyantoro, Handayani 2017). In other words, most forests in the coastal dune are human-made and thus included in the cultural, instead of biotic, element. An example of a biotic element is Borassus flabellifer, which is a high, robust tree that is sparsely distributed and the covered area does not meet the minimum mapping unit.

The research has added another class to describe the human interventions in the study area, that is, info-telecommunication. It consists of transceiver and monitoring stations whose area coverage meets the minimum mapping unit. These constructions are analysed in LER because their height can create barriers to sand-carrying wind and influence the coastal dune formation. However, the structures generally have voids of varying sizes, creating a surface through which wind can still flow, and for this reason, they are not included in the urbanogenic landscape.

The landscape ecological risk (LER) values, as expressed in ERI, range from 0.0 to 1.0. This range of values likely results from the selected spatial interpolation method or the multiplication of landscape disturbance and fragility. In the spatial interpretation, the sampling points are grids of equal size with a centroid, resulting in an equidistant ERI. Kriging indicators are believed to contribute to the derived index values (Kumar, Sinha 2018). Because the sampling points are distributed evenly, the spatial value of the LER generates a very low RMSE. The four scenarios have a somewhat wide variety of landscape disturbance index (0.37–68.18) and landscape fragility index (0.048–1.0). Because the disturbance and fragility index values were multiplied proportionally within the grid, the LER value (or ERI) can vary from 0.0 to 1.0.

Based on the derived LISA cluster maps, bare land, grassland and sparse shrubs form high-high clusters. These three sub-classes naturally have a high landscape fragility index that contributes to the LER calculation (Zhang et al. 2018). Sheds, non-irrigated croplands, roofed/ indoor tourist attractions and parking areas also form high–high clusters, which is associated with high landscape disturbance (Liu et al. 2020), especially intensive fragmentation (Wang et al. 2020). Meanwhile, low–low clusters can be found in rice fields, non-irrigated croplands, livestock farms and sparse shrubs with high fragility but low fragmentation degrees. In addition, dense shrubs, multiple-species forests, settlements, Casuarina forests, roofed/indoor tourist attractions and unroofed/outdoor tourist attractions are LULCs with low–low clusters because of the relatively low fragility.

The derived LISA cluster maps also show that low–low clusters are mainly found in natural, tree agrogenic and urbanogenic landscapes. Tree agrogenic and urbanogenic landscapes are slightly fragile because of the low possibility of being converted to other landscapes, whereas the natural landscape is highly fragile. Previous studies assume that fragmentation contributes more to landscape disturbance than dominance and isolation (Xie et al. 2013, Zhang et al. 2013, Liu et al. 2020). Because the natural landscape is the matrix and the largest patch of the coastal dune observed, its low LER value is attributed to the low degree of fragmentation (Li et al. 2020). Even though it is slightly fragile, the tree agrogenic landscape also forms high–high clusters because of the high disturbance in the northern supporting zone. Other than the tree agrogenic landscape, high–high clusters are distributed in agrogenic, industrogenic and tourism-sports landscapes chiefly because these types of intervention are relatively highly fragile.

There are some differences in LER between the class and sub-class scenarios. The sub-class scenario has bare land as the matrix or the largest patch, while the class scenario has natural landscape. Bare land is a sub-class included in the natural landscape, which has the highest landscape fragility (Zhang et al. 2018). On the one hand, the bare land forms a high–high cluster; on the other hand, the natural landscape is a low–low cluster. Despite the same level of fragility, their disturbance levels are different because of the landscape fragmentation: high in bare land (this sub-class has the greatest number of patches) and low in natural landscape because of the merging of bare land, grassland, sparse shrubs, litter accumulation, footpaths and natural river. Other than these differences, the LER levels derived from the class and sub-class scenarios generally form a similar pattern, i.e., none of the scenarios indicate a particular risk level that represents a linear landscape structure.

Figure 8 can also be used to explain Figure 7 to obtain information about the temporal variation of LER clusters. There are more LER clusters in 2020 than in 2011 for both class and sub-class scenarios. However, based on the size, the 2020 clusters are generally smaller than the 2011 clusters. An increase in the number of clusters is in line with the addition of patches, whereas a decrease in the cluster size indicates a reduction in the mean size of the patch. According to Forman (1995), patch number addition and mean patch size reduction typify fragmentation and dissection (fragmentation by linear elements). Because the LER levels show no indication of risk from a linear landscape structure, it can be inferred that the spatial ecological process occurring in the study area is fragmentation.

Land transformation caused by human interventions is associated with the interaction of economy and ecology. Although human intervention tends to prioritise economic interests, it is not always destructive but can give some positive values (Kubalikova et al. 2019). Landscape has multiple functions including ecological, social and economic (Li et al. 2020). Therefore, it is imperative that spatial planning to be applied (Cao et al. 2019, Wang et al. 2020) considering factors in natural habitats and ecosystems (e.g., bare land or natural landscape) as primary conservation areas (Liu et al. 2020). The Parangtritis coastal dune has been restored to protect the bare land or natural landscape (Hendrastuti et al. 2018), even though the restoration initiatives often focus on arranging fisheries and built-up areas (Laily et al. 2018). Restoration was selected as the best way to manage the Parangtritis coastal dune. It will help the sand dunes back to its initial condition slowly by reducing the remaining trees planted in the program since 1975 (Fakhruddin et al. 2010), without any hard structures. It does not prohibit people to do restricted economic activity. Some non-native flora and fauna are still alive on it (Łabuz 2015).

For ecological protection, it has been proposed to prohibit afforestation and construction in the Parangtritis coastal dune, especially in the wind tunnel (Sunarto 2014), which includes the foredune and interdune of the core zone. However, in some other studies, reducing ecological risk means improving the environment through reforestation (Peng et al. 2015, Zhang et al. 2018, Gong et al. 2021) because high LER is frequently assumed to result from deforestation (Wang et al. 2020) and desertification (Xie et al. 2013). Furthermore, strategies to reduce LER are optimal when aimed at increasing matrix connectivity (Xie et al. 2013).

This study has several limitations. It uses SFAP as reliable spatial data in a fine resolution to accommodate detailed mapping and LEU classification. However, for individual research, SFAP can be costly. Moreover, the most recent technology such as artificial intelligence has not been applied together with SFAP, which can otherwise facilitate the delineation of LEU and substantially reduce the time spent to perform visual interpretation. Human intervention is complex, and some may not be observable by SFAP. Furthermore, the study uses a relatively simple calculation method (Cao et al. 2019), resulting in LER values that are sensitive to the grid size (Zhang et al. 2018) and scale (Liu et al. 2020). The LER assessment is also to some extent subjective; therefore, AHP was used in the landscape fragility assessment to reduce the subjectivity. The derived model is static (Cao et al. 2019, Li et al. 2020) and only analyses the available data to illustrate the current condition, which is contrary to disaster risk assessment that projects future possibilities with existing data. The LER assessment results are also uncertain (Cao et al. 2019; Liu et al. 2020) because the evidence in the field is rather lacking and is only limited to LULC and manifestation of human intervention. To provide recommendations for spatial planning, the LER assessment needs to use more reliable and specific methods to produce distinct results.