GIS-based landform classification of Bronze Age archaeological sites on Crete Island

Various physical attributes of the Earth’s surface are factors that influence local topography and indirectly influence human behaviour in terms of habitation locations. The determination of geomorphological setting plays an important role in archaeological landscape research. Several landform types can be distinguished by characteristic geomorphic attributes that portray the landscape surrounding a settlement and influence its ability to sustain a population. Geomorphometric landform information, derived from digital elevation models (DEMs), such as the ASTER Global DEM, can provide useful insights into the processes shaping landscapes. This work examines the influence of landform classification on the settlement locations of Bronze Age (Minoan) Crete, focusing on the districts of Phaistos, Kavousi and Vrokastro. The landform classification was based on the topographic position index (TPI) and deviation from mean elevation (DEV) analysis to highlight slope steepness of various landform classes, characterizing the surrounding landscape environment of the settlements locations. The outcomes indicate no interrelationship between the settlement locations and topography during the Early Minoan period, but a significant interrelationship exists during the later Minoan periods with the presence of more organised societies. The landform classification can provide insights into factors favouring human habitation and can contribute to archaeological predictive modelling.


Introduction
During the Middle Minoan period of the Bronze Age period in Crete, there appears to have been a widespread increase of sites in low elevation areas suitable for arable farming. Systematic archaeological surveys across eastern and central Crete have examined the settlement dynamics and in all cases highlight a generalised population movement from high elevation areas with limited arable land to lower elevation areas, particularly plains favourable to cultivation and efficiency in irrigation [1]. Most of the studies hypothesized that this population movement tendency was caused by economic and political reasons but only a few of them considered the possibility that this tendency could be a result of environmental conditions and a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 complex rough terrain, mainly hills and ridges [6]. The Mirabello Bay area consists of small coastal valleys: most of the settlements seem to be placed on slopes rather than on the best (flat, lowland) agricultural land [42]. The site pattern was defensive, providing lines of sight between the settlements and out across the Mirabello gulf. Thus the Vrokastro site pattern may be in some accord with [43] suggestion that "people were forced to look for their safety in barely accessible places". During this period it seems there was a dichotomy between occupying lowlying coastal sites and upland locations, indicating a tendency of the population to occupy diverse types of terrain, seeking to gain control of water sources and sites with the greatest visibility [6]. The evidence for settlement hierarchy during this period has been doubted [44], with clusters of settlements found around the arable coastal basins and not around large settlements.
The Middle Minoan is characterized by settlements in homogeneous landscapes, such as the Kavousi plains, which were well-suited for cultivation [38]. There appears to have been a movement of population into the lowlands and arable upland areas close to water resources, with an associated increase in number of settlements near the best agricultural land during Middle and Late Minoan [6,33,38]. [5] notes that the pattern of settlements exploitation in the Kavousi area is in accordance with the Argolid model during this period, with concentrated population around arable areas, similar to the expansion pattern of small farms inland in Chania district.
Between the Early and Middle Minoan, population growth occurred in the Gournia valley, with continuous expansion of settlements and coastal trading interests [6]. In addition, during the Middle Minoan there is exploitation of the landscape at relatively high elevations, between 450 m and 700 m above sea level. River terraces and bedrock of conglomerates and marl provided fertile soil for cultivation of the Gournia valley, along with fault fractures that provided water access via springs at higher elevation [6].
In the Late Minoan an overall reduction of the number of settlements is observed, especially in upland sites. Despite that, rural settlement is oriented to arable land and water resources, preferentially between 50-200 m elevation range, while a retreat from coastal zones (0-50m) seems to have occurred [6]. An event that might be able to explain this reduction of settlements is the impact of the Theran eruption, with evidence of a thick layer of tephra/ash found in the Mochlos area, near Gournia. Population centralization is also observed during the Late Minoan, for example at Pyrgos or Gournia. Some new large building structures during this period can be found (e.g. around the Istron river valley) which appear to be linked to trade routes, access to arable land and water resources [45].

Materials and methods TPI and DEV analyses for various neighbourhood sizes
The landform classification was based on TPI and DEV analyses, using as an input dataset the free ASTER Global DEM (G-DEM: 30m pixel size) and the allocated archaeological sites for each Minoan period, as recorded by the archaeological surveys of [6], [33] and [38] (Figs 2 &  3). A buffer zone of 300m around each site was selected for the analyses, in order to evaluate the surrounding dominant landscape of each individual site. The TPI and DEV analyses were evaluated for various neighbourhood sizes, to select the optimum neighbourhood for the rest of the methodological framework. Such a decision depends on the degree that a particular neighbourhood size will represent, or separate better the various morphological classes in a more coherent and uniform way. The Fragstats free software was used to calculate the patch density (PD) and numbers of patches (NP), in order to check the fragmentation of the patches for the various morphological classes ( Table 1). The NP equals to the number of patches of the  GIS-based landform classification of Bronze Age archaeological sites, Crete corresponding patch type, as a measure of the extent of fragmentation of the patch type; and PD expresses in addition the number of patches on a per unit area basis [46].
TPI or DIFF (z o -z) measures the relative topographic position of the central point as the difference between the elevation of this central point and the mean elevation within a pre-determined neighbourhood [30,31]. The result from TPI is a positive value when the central point is situated higher than its average neighbourhood, while a negative value indicates a lower location than the average neighbourhood. This index can be used to classify the landscape into morphological classes [47] (Fig 2).
DEV ( zoÀ z SD ) measures the relative topographic position of the central point as the TPI divided by the standard deviation of elevation (SD), within a predetermined neighbourhood, where SD measures the variability of the cell values in a DEM, around this central point, within the

Indices
Description Range

SHAPE
Equals patch perimeter (m) divided by the square root of patch area (m 2 ), adjusted by a constant to adjust for a square standard.

SHAPE≦1,without limit.
SHAPE = 1 when the patch is square and increases without limit as patch shape becomes more irregular.

PROX
Equals the sum of patch area (m 2 ) divided by the nearest edge-toedge distance squared (m 2 ) between the patch and the focal patch of all patches of the corresponding patch type whose edges are within a specified distance (m) of the focal patch. When the search buffer extends beyond the landscape boundary, only patches contained within the landscape are considered in the computations. Note that the edge-to-edge distances are from cell center to cell center.
PROX!0, PROX = 0, if a patch has no neighbors of the same patch type within the specified search radius. PROX increases as the neighborhood (defined by the specified search radius) is increasingly occupied by patches of the same type and as those patches become closer and more contiguous (or less fragmented) in distribution. The upper limit of PROX is affected by the search radius and the minimum distance between patches.

LSI
Equals .25 (adjustment for raster format) times the sum of the entire landscape boundary (regardless of whether it represents 'true' edge or not, or how the user specifies how to handle boundary/ background) and all edge segments (m) within the landscape boundary involving the corresponding patch type, divided by the square root of the total landscape area (m 2 ). Note, total landscape area (A) includes any internal background present.

LSI!1,without limit.
LSI = 1, when landscape consists of a single square patch of the corresponding type and increases without limit as landscape shape becomes more irregular.

CIRCLE
Equals 1 minus patch area (m 2 ) divided by the area (m 2 ) of the smallest circumscribing circle.

0<CIRCLE<1,
CIRCLE approaches 0 for circular patches and approaches 1 for elongated linear patches.
COHESION Equals 1 minus the sum of patch perimeter (in terms of number of cell surfaces) divided by the sum of patch perimeter times the square root of patch area (in terms of number of cells) for patches of the corresponding patch type, divided by 1 minus 1 over the square root of the total number of cells in the landscape, multiplied by 100 to convert to a percentage. Note, total landscape area (Z) excludes any internal background present.

0<COHESION<100,
COHESION approaches 0 as the proportion of the landscape comprised of the focal class decreases and becomes increasingly subdivided and less physically connected. COHESION increases as the proportion of the landscape comprised of the focal class increases until an asymptote is reached near the percolation threshold.
PAFRAC Equals 2 divided by the slope of regression line obtained by regressing the logarithm of patch area (m 2 ) against the logarithm of patch perimeter (m). That is, 2 divided by the coefficient b1 derived from a least squares regression fit to the following equation: ln(area) = b 0 + b 1 ln(perim). Note, PAFRAC excludes any background patches.
A fractal dimension greater than 1 for a 2-dimensional landscape mosaic indicates a departure from a Euclidean geometry (i.e., an increase in patch shape complexity). It approaches 1 for shapes with very simple perimeters such as squares, and approaches 2 for shapes with highly convoluted, plane-filling perimeters.

PD
Equals the number of patches of the corresponding patch type divided by total landscape area (m 2 ) PD>0, constrained by cell size.
PD is constrained by the grain size of the raster image, because the maximum PD is attained when every cell is a separate patch. Ultimately cell size will determine the maximum number of patches per unit area predetermined neighbourhood [30]. The results from DEV are: positive, when the central point is situated higher than its average neighbourhood; or negative, when the central point is situated lower than its average neighbourhood (Fig 3). The exact difference in height between the central point and the average height of the neighbourhood can be highlighted by those analyses, with the range of the outcome values being dependant on the elevation differences in the DEM within the selected neighbourhood. The selection of the neighbourhood sizes were based on 5 out of 12 candidate radii from 100 to 2000m and the ones used to [19] study were decided as optimum for this study. The circular neighbourhood sizes considered within this study as more useful had radii of: 100 m, 300 m, 600 m, 1200 m and 2000 m. The application of a neighbourhood size is important during the analysis and is related to the landscape feature being analysed. In order to classify small features such as streams and valleys a small neighbourhood size need to be used. To identify large canyons or other large landscape features a large neighbourhood size should be selected during the analysis. When we are dealing with larger size of the neighbourhood and large differences in elevation values these conditions result to the wider range of the output values. Based on the above observations TPI was used instead of DEV for the subsequent analyses of slope position and landform classification, as it represents a more regional representation of the morphological classes (

Slope position classification for various neighbourhood sizes
In order to classify the landscape into slope position classes, the variation of TPI values by the slope at each point was examined. High TPI values determine the tops of the hills, while low TPI values highlight valley bottoms. TPI values close to 0 are observed over flat to mid-slope areas. In this study, a six-category slope position class ( Table 2) was determined for each of the five neighbourhood sizes considered (100 m, 300 m, 600 m, 1200 m and 2000 m) (Fig 4). A 5s lope threshold value was considered for the discrimination of flat areas and mid-slope areas, while a TPI threshold value of ±1SD was considered for the identification of hilltops and valley bottoms. The six slope position classes, based on [31], are presented in Table 2.

Landform classification for various neighbourhood sizes
The landform classification was based on TPI analysis. The TPI based landform classification can produce 10 landform classes: streams, mid-slope drainages, local ridges, valleys, plains, foot slopes, upper slopes, upland drainage, mid-slope ridges and high ridges [21] (Fig 5). Usually, two neighbourhood sizes are combined to offer detailed geomorphological information through the discrimination of complex landscape features, as a single neighbourhood size provides less information about the general shape of the landscape [21]. In this study, the two combined neighbourhood sizes in each case were: i) 100 m and 600 m; ii) 300 m and 1000 m; iii) 300 m and 2000 m and; iv) 600 m and 2000 m. Table 2. Classification of the landscape into morphological classes, where SD is Standard Deviation.

Morphological classes Weiss (2001) [31]
Valley bottoms TPI<-SD The allocated settlements for each Minoan period in Vrokastro district provided an accuracy assessment of the combined neighbourhood sizes outcomes. Vrokastro district was selected to test the methodology on a more heterogeneous landscape relative to the less rough terrain of Phaistos region. The accuracy assessment consisted of the comparison between the coverage percentage of the individual landforms classes based on TPI analysis and the geological-geomorphological description of the settlements provided by the archaeological surveys of [6]. The Fragstats freeware was used to determine which of the combined neighbourhood sizes provided the best representation of the various landform types, relative to the rest of the tested neighbourhood sizes, regarding their shape and fragmentation characteristics. Various indices were examined via Fragstats (Table 1): i) shape index (SHAPE); ii) proximity index (PROX); iii) landscape shape index (LSI); iv) related circumscribing circle (CIRCLE); v) patch cohesion index (COHESION); vi) patch density (PD); vii) numbers of patches (NP) and; viii) perimeterarea fractal dimension (PAFRAC).
A spatial autocorrelation evaluation based on Moran's I approach was carried out for the EM, MM, LM settlements and 100 random sites of the study area, to check whether the pattern expressed was clustered or random [48]. The derived z-score or p-value can indicate statistical significance, with a positive Moran's I index value indicating tendency toward clustering, while a negative Moran's I index value indicates tendency toward dispersion. The z-score and p-value calculations indicate whether the null hypothesis can be rejected. In this case, the null hypothesis states that features are randomly distributed across the study area. Clustering is a statistical classification technique for separating spatial information into relatively homogeneous groups: similarities are high between members belonging to a class, or cluster [48].

TPI and DEV analyses for various neighbourhood sizes
The TPI and DEV outcomes are presented on Figs 2 and 3 for the various neighbourhood sizes (150 m, 300 m, 600 m, 1200 m and 2000 m), with six slope classes being discriminated ( Table 3). The results show that TPI provides a generalised representation of topography, while DEV highlights subtle topographic variations, which is in accordance with the findings of [32]. This study focuses on the regional representation of the morphological classes so the generalised approach of TPI was selected, instead of DEV, for the analyses of slope position and landform classification evaluation. The calculated indices of NP and PD, revealed that DEV provides a more fragmented representation of the morphological classes than TPI, with the PD values indicating a higher density of patches for the individual classes ( Table 3). The comparison of the NP and PD values between the various neighbourhood sizes, for both DEV and TPI, indicated that the 150 m neighbourhood size shows high fragmentation with a large number of patches and high density of patches. On the other hand, the 600 m and 1200 m neighbourhood sizes highlight a uniform composition of the morphological classes (steady decrease of NP and PD), which results to a high degree of generalization regarding the morphological classes, with particular features not visible at that scale, due to the generalization ( Table 3). As a result, the 300m neighbourhood size has been selected as optimum for the rest of our analyses: it discriminates the various features with less fragmentation and diverse NP and PD values, without a high degree of generalization.

Slope position classification for various neighbourhood sizes
Slope position classification based on TPI for Early, Middle and Late Minoan periods consisted of six classes (Fig 4). The percentages of the slope position classification areal coverage, for the individual Minoan periods, are presented on Table 4. For valley bottoms, a decrease in the coverage percentage exists with increasing neighbourhood sizes, for all Minoan periods. For lower slopes, there is an increase of coverage with a neighbourhood size increase, for all Minoan periods. Gentle slopes, steep slopes and ridges have a steady variable coverage percentage for all Minoan periods. The two notable observations on those landform types are that: i) the upper slopes reveal a high decrease between 100 m and 300 m neighbourhood sizes, for all Minoan periods, but then a steady variation exists; ii) during the Early Minoan period, the ridges class coverage increases gradually with increased neighbourhood sizes; while during the Middle and Late Minoan periods there is a steady variation of ridges class coverage between the neighbourhood sizes. This might be a result that during Early Minoan period the settlements were developed on high hilltops, presumably for defensive purposes, so they were located on rougher terrain relative to Middle and Late Minoan. The neighbourhood size of 600m provides the better option to characterize the region, based on the indices statistics, visual interpretation of satellite imagery, in conjunction with information from the derived geomorphometric variables and the landscape shape/fragmentation indices (NP, PD, PROX and LSI) (Fig 4; Tables 4 & 5). The NP and PD indices indicate that lower neighbourhood size represents a fragmented representation of the morphological classes, with a higher density of patches for the individual classes. As the neighbourhood size increases, a uniform composition of the morphological classes (steady decrease of NP and PD) is observed. These indices, in conjunction with the rest of the indices, imply that the 600 m neighbourhood size is the optimum neighbourhood (Table 5). This particular neighbourhood size defines better the morphological classes, giving equal emphasis on regional and subtle topography, without exaggerating one or other morphological classes as it happens with the rest of the neighbourhood sizes.

Landforms classification for various neighbourhood sizes
The TPI based landform classification for the Early, Middle and Late Minoan periods consisted of ten different landform types, for combined neighbourhood sizes (Fig 5). The combined neighbourhood sizes of 300 m and 1000 m were optimal for characterizing the region, based on the Fragstats indices statistics and visual interpretation of satellite imagery (Table 6). Moreover, the study areas are located on rough terrain with abrupt changes of the topography, so a more generalised overview from the combined neighbourhood sizes of 300 m and 1000 m is considered as optimum for larger scale geo-archaeological interpretations. The more subtle topography determined by the combined neighbourhood sizes of 100 m and 600 m would be suitable for flat areas/benches, with the landform types being less fragmented. The histogram of landform elements percentage for the combined neighbourhood sizes of 300 m and 1000 m reveals the predominant landform types during each of the Minoan periods (Fig 6).
For the Kavousi-Vrokastro district, two main highlights of the settlements location, from the Early Minoan towards the Late Minoan are: i) high increase on deeply incised streams, valleys and plains; ii) high decrease on open slopes, mid-slope drainages and mid-slope ridges (Fig 6). The settlement types and the number of settlements are presented in Table 7. For the Phaistos district, two main highlights of the settlements location from the Early to the Late Minoan are: i) decrease of high ridge and upper slope locations; ii) steady increase of locations in deeply incised valleys, plains and open slopes (Fig 6). The Phaistos district has less rough terrain than Kavousi-Vrokastro, so the tendency of populations moving to flatter landscape from the Early to Late Minoan, is more pronounced. On the other hand, there is a drop in the number of settlements found on high ridges, which indicates a tendency towards abandoning high-elevation settlements. These interpretations are in accordance with findings from the archaeological surveys of [6,33,38]. Based on Moran's I analysis of the EM, MM, LM settlements and 100 random sites to check whether the pattern expressed is clustered or random, the outcomes from the Kavousi-Vrokastro study showed the following case study: i) for the 100 random points the z-score was -0,75 and p-value 0.45, indicating that the pattern is random; ii) similar observations exist for the EM settlements with the z-score being 1.05 and p-  value 0.29. However, for the MM settlements the z-score was 7,76 and p-value 0, indicating that there is less than 1% likehood that this clustered pattern could be a result of random chance. Furthermore, with LM settlements the z-score was 3.15 and p-value 0,001, indicating a clustered pattern. Similar observations also exist for the Phaistos case study: i) for the 100 random points the z-score was 1.47 and the p-value was 0.14, indicating that the pattern does not appear to be significantly different than random; ii) similarly random pattern exists for the EM settlements with the z-score being -0.98 and the p-value being 0.32. On the other hand, for the MM settlements the z-score was 5 and the p-value was 0, indicating that there is less than 1% likehood that this clustered pattern could be a result of random chance; while for the LM settlements the z-score was 2.84 and the p-value was 0,0044, indicating a clustered pattern.

Discussion
This study has quantified the spatio-temporal variations of Early, Middle and Late Minoan settlement locations in the landscapes of Crete. A general trend is observed during the Early to Middle Minoan, with the population moving inland in search of arable land, most settlements initially being located within 1.5 km of the coast. During the Middle Minoan, there was an increase of settlements over low-gradient slopes, upper slopes and mesas in hilly terrain, where they remained during the Late Minoan, probably because of their strong defensive advantage (Fig 6). The analysis of settlement areas indicates that during the Middle Minoan there was: i) an increase of settlements in deeply incised valleys, with people looking for access to perennial water supplies and; ii) a reduction of the number of settlements in shallow U-shaped valleys, mild slopes & mid-slope ridges, perhaps due to movement of population to more arable land on the plains and U-shaped valley floors (Fig 6) [6,33,38].
In Phaistos district, the Early Minoan settlements are located in heterogeneous landscapes, distributed sparsely over variable terrain. In the Middle Minoan the settlements are found over homogeneous landscapes, such as the plains. In the Late Minoan the largest settlements are found on the plains and lowlands. During the Middle to Late Minoan, the number of settlements decreased, especially at higher elevations; the settlements on the plains remained, while farming increased and population concentrated in the larger settlements [6,33,38].Based on the analysis of the area percentage of the Minoan settlements in Phaistos district, population increased within the deeply incised valleys, perhaps due to the need for perennial water, and on the plains and foot slopes, perhaps due to increased arable land access. However, population decreased on the upper slopes, mesas and high ridges, perhaps as a result of the observed movement to arable lowlands (Fig 6).
The accuracy assessment, comparing the findings of this study with the archaeological surveys of [6] (Table 8), shows a low percentage of no agreement (14%) while there is a substantial

Ioannimiti region
These sites are located on the hilltops, slopes and ridges. IM1 located on an east facing rise, IM2 lies along the western side of the Vathi peak, with excellent view of Ioannimiti area. Terrace walls constructed of large boulders still stand on the gentle, 15 0 , south facing slopes where IM3 is located. IM9 on a low south facing slope above a torrent bed at the centre of the promontory.

Conclusion
During the last decade the quantification of landform types has attracted the interest of the geoinformatic research community, with various studies using GIS-based approaches. The evaluation of geomorphometric datasets can be integrated with GIS techniques, highlighting the information within the interlinked geographic data. The extracted information can be particularly useful for landform classification, as this case study of Crete has illustrated. Such information is useful when linked to geospatial data about the distributions of archaeological sites over space and time. This investigation of Phaistos, Kavousi and Vrokastro districts has produced valuable information regarding the distribution of settlements during the Early, Middle and Late Minoan periods. Based on the analysis of this study, a general trend is observed during the Early to Middle Minoan, with population moving from heterogeneous terrain at higher elevation, to the lowlands (Fig 6). During the Middle to Late Minoan, the population remained in the arable lowlands, with better organization and concentration in larger settlements.
This study was constrained by the limited amount of paleo-environmental data available for Crete, resulting in a knowledge gap that adds to the uncertainty associated with the interpretations made about the factors driving the variations in the Minoan settlement distributions. Nevertheless, the methodology presented here can provide useful spatio-temporal analyses, at district scales, for future studies to examine at local scales, with associated studies of palaeo-environmental conditions or archaeological predictive modelling. In addition, this study offers valuable information for further research, where socio-economic or political factors can be considered for settlement hierarchy assessments.