High and Medium Resolution Satellite Imagery to Evaluate Late Holocene Human–Environment Interactions in Arid Lands: A Case Study from the Central Sahara

We present preliminary results of an Earth observation approach for the study of past human occupation and landscape reconstruction in the Central Sahara. This region includes a variety of geomorphological features such as palaeo-oases, dried river beds, alluvial fans and upland plateaux whose geomorphological characteristics, in combination with climate changes, have influenced patterns of human dispersal and sociocultural activities during the late Holocene. In this paper, we discuss the use of mediumand high-resolution remotely sensed data for the mapping of anthropogenic features and paleoand contemporary hydrology and vegetation. In the absence of field inspection in this inaccessible region, we use different remote sensing methods to first identify and classify archaeological features, and then explore the geomorphological factors that might have influenced their spatial distribution.


Introduction
Analysing spatio-temporal characteristics of human-driven landscape change is key for assessing the ecological and socioeconomic consequences of landscape transformations.Deserts have been-and still are-places of innovation, where humans have been fluidly adapting to extreme environmental conditions.Identifying evidence of human activity in arid lands is crucial to understand how past (and present) communities coped with the scarcity of natural resources.Accessibility of many North African and Middle East desert regions has always been limited, due to their physiography and climate.In recent years, access has been further reduced due to political instability, civil war and unsafe fieldwork conditions.As a consequence, archaeological investigations in these regions have been profoundly affected, with the majority of field-based and long-term research projects interrupted.In general, archaeological research in North Africa has focused on the so-called "Green Sahara", i.e., the wetter period from ca. 15,000 to 5000 years before present (BP), when major socioeconomic transformation occurred in the region, including the shift to food production through the adoption of pastoralism, and the early use of cattle by-products.Late prehistoric, historical, and even modern and current activities, which correspond to the last ~5000 years (including the time periods of the "Dry Sahara" or "Brown Sahara"), are less well-known.Collecting information from cultural landscapes that are in danger of disappearing due to mining and oil/gas drilling [1][2][3], agricultural expansion [4,5] and, more recently, conflict [6] is paramount, not only for the preservation of the historic memory of these landscapes and in disentangling their dynamics, but also for future planning, conservation and development.
Remote sensing (RS) and Geographic Information Systems (GIS) have demonstrated great potential for creating spatio-temporal datasets at a relatively low cost and high speed, as compared to field data collection [7,8].This is particularly true in hyper-arid Africa, where collected field data are scarce, the areas to be covered are large, and the visibility of certain classes of archaeological evidence is high, since it contrasts with the background barren environment.In order to refine high temporal and spatial resolution datasets that may serve as background for the study of human-environment interactions in the Sahara, we have developed a combined methodology for mapping a number of indicators of, and drivers for, human occupation in this environment, in particular during the "Dry Sahara" period.This procedure allows for relatively rapid regional-scale archaeological mapping of visible remains, but also identifies proxies of ancient and current environmental settings where human occupation may have been concentrated in the past, such as river networks (paleohydrological mapping) or in humid areas (based on vegetation mapping of springs or oases).This procedure is aimed at creating an analytical environment where spatial analysis and spatial modelling approaches can be applied to understand the evolution of desert lifeways during the late Holocene.
The Central Sahara: Archaeology, Environment, and Remote Sensing In the Central Sahara, previous work on the archaeology of the last 5000 years has focused mainly on SW Libya.The archaeological evidence in this region is represented by (i) rock shelters and caves, and (ii) open air sites.The records of rock shelters and caves have played a crucial role in the reconstruction of Holocene cultural trajectories [9], because their deposits host long sequences of Holocene occupation, from Epipalaeolithic hunter-gatherers to Neolithic pastoralists.However, in many cases the upper strata of archaeological deposits in rock shelters have been removed by wind.A few exceptions are some rock shelter sites in the Libyan Tadrart Acacus that have yielded 14 C dates from occupation horizons of ca.3000-1000 years BP [10].On the other hand, recent studies of rock art [11] and carved texts [12,13] from SW Libya have demonstrated how a landscape approach to this evidence can shed new light on late prehistoric, historical, and modern human-landscape interactions.Most archaeological data, however, come from open-air sites that have higher visibility in the landscape.These sites include temporary pastoral campsites, stone cairns, villages, and outposts.Different types of Saharan funerary monuments have been described [14][15][16][17][18][19] and some new radiometric dates are available from recent excavations.Study of Garamantian (ca.1000BC-AD700) forts, outposts, and compounds [4,5,[20][21][22] has shed light on a complex network of trade within the region, that further developed in medieval times [23,24].In spite of this body of work, knowledge of the distribution of archaeological sites from the last 5000 years in the central Sahara is patchy, although recent work has started systematically mapping and cataloguing sites that are visible from high-resolution satellite imagery [4,5,25].Despite the lack of refined chronological assessment and inevitable de-contextualization of the archaeological record, this approach has populated large "empty" spaces, enabling systematic analysis of site categories and spatial distributions and setting a powerful basis for future field and remote investigations.
The Central Sahara region includes a variety of active and relict geomorphological features, such as palaeo-oases, dried river beds, lowland alluvial fans, and upland plateaux whose geomorphological characteristics, in combination with climate changes, have influenced patterns of human dispersal and socio-cultural activities during the Holocene (e.g., [26]).Landscape characterization and geomorphic mapping are therefore fundamental for decoding the archaeological data of various types.In this paper, we discuss the use of different mid-and high-resolution remotely sensed data for the mapping of archaeological evidence, and landscape palaeo-and contemporary hydrology, and vegetation.
Palaeohydrology has proven to be key in understanding the dispersal of hominins in Africa and the Arabian Peninsula during the Quaternary, and establishing where and why humans have moved in and out of arid areas [27,28].Extensive and interconnected hydrological systems facilitated human dispersal during the "Green Sahara" phase [27] and also played an important role during drier periods when wadi beds were the focus of human activities and movements, as they are under similar conditions today [29].The identification of river valleys using digital topographic data is therefore the first step to characterize Saharan paleohydrology, and the production of a high-resolution digital drainage network has been one of the first steps in this project.
Digital hydrological networks such as HydroSHEDS, a global hydrologically conditioned dataset based on Shuttle Radar Topography Mission (SRTM) data, are publically available [30].Nevertheless, such global products are released at a scale that fails to capture the smaller channels that are necessary for a sufficiently detailed representation of the region's hydrology.This limitation can be overcome by using multiple digital elevation models (DEMs) and through standard digital hydrological network extraction techniques included in most off-the-shelf GIS software packages.
In the Central Sahara, riverbeds host subsurface and surface moisture and are the foci of past and present human occupations.Given there is little difference in the climatic conditions between the full establishment of aridity in the region 5000 years ago and today [31], the mapping of the current vegetation pattern can be used to infer past conditions, and in turn indicate those areas that would have remained habitable during previous dry periods.Contrary to the popular belief that considers desert oases as the only areas with vegetation, the Central Sahara hosts a variable but permanent plant cover, along with an ephemeral grass cover that may be present across large regions following rain [32,33].Drought-tolerant trees and shrubs represent an important natural resource, due to the high degree of water retention in their plant tissues.This is key to current pastoral systems in Saharan mountains, as these plants provide nourishment to sheep and goats.The trees and shrubs are also a clear indicator of subsurface humidity, and can be used as a guide to identify suitable hotspots of past human activity.
It is widely accepted that remote sensing is a significant data source to study Land Use/Land Cover (LULC) and vegetation changes [34] in areas with high precipitation variability.Two main approaches have been used for this purpose: vegetation indices and classifications.Vegetation indices such as the Normalized Difference Vegetation Index (NDVI) exploit the high reflectance of vegetation in the near-infrared and high absorption in the red bands of the spectrum.The contrast between these channels can be used as an indicator of vegetation health, since it is a parameter that correlates well with photosynthetic activity [35].NDVI is also a good proxy for the presence of vegetation in an area, since NDVI values between 0.3 and 0.8 are an indication of vegetation, whereas values of 0.2 to 0.3 represent bare soils, and negative NDVI values represent water.Feature extraction from multispectral remote sensing images is thus possible by thresholding NDVI values, and can be used spatially to map different LULC classes in a given area.The NDVI index has also been used to detect high water table zones in humid warm-temperate regions [36] and to evaluate the effects of precipitation in semi-arid landscapes [37].
Classification is the process of extracting differentiated classes or themes (such as land use categories or vegetation species) from raw remotely-sensed satellite data.Unsupervised classification has been traditionally used for thematic mapping of areas for which no ground data are available, while supervised classification requires training datasets, generally collected through ground surveys, which are then used as predictor variables to assign pre-determined known classes to the sampled units.Maximum likelihood classifiers (MLC) have been extensively used for supervised classification methods but, since they require the spectral response of each class to follow a Gaussian distribution, their results are limited in quality, due to the fact that most data (in particular in complex landscapes) do not follow such distributions, meaning that the accuracy of class separation is reduced [38].Non-parametric classifiers, that do not rely on normal distributions, such as support vector machine and random forest classifiers, have been shown to outperform MLC in many settings [39][40][41].
In this paper, we present the methodology and first results of an integrated study of the human (anthropogenic) and geomorphic landscape elements from the Central Sahara, mapped using remote sensing data and analysed using GIS approaches (Figure 1).This paper focuses exclusively on the remotely-sensed proxies that can be used to develop a more nuanced understanding of human-landscape interactions during the last 5000 years in the region.Satellite data analyses at varying resolutions coupled with spatial analyses of topographic derivatives and environmental proxies using GIS provide a powerful approach to the mapping of large areas in a semi-automated manner.In this paper, we introduce the study area, then outline the datasets and methods used (for land cover, palaeo-hydrology, and identifying anthropogenic signatures), and then provide some local examples illustrating these human-landscape relationships, based on detailed numerical analysis and modeling from remote sensing data coupled with contextual archaeological data.
Remote Sens. 2017, 9, 351 4 of 21 are available, while supervised classification requires training datasets, generally collected through ground surveys, which are then used as predictor variables to assign pre-determined known classes to the sampled units.Maximum likelihood classifiers (MLC) have been extensively used for supervised classification methods but, since they require the spectral response of each class to a Gaussian distribution, their results are limited in quality, due to the fact that most data (in particular in complex landscapes) do not follow such distributions, meaning that the accuracy of class separation is reduced [38].Non-parametric classifiers, that do not rely on normal distributions, such as support vector machine and random forest classifiers, have been shown to outperform MLC in many settings [39][40][41].
In this paper, we present the methodology and first results of an integrated study of the human (anthropogenic) and geomorphic landscape elements from the Central Sahara, mapped using remote sensing data and analysed using GIS approaches (Figure 1).This paper focuses exclusively on the remotely-sensed proxies that can be used to develop a more nuanced understanding of human-landscape interactions during the last 5000 years in the region.Satellite data analyses at varying resolutions coupled with spatial analyses of topographic derivatives and environmental proxies using GIS provide a powerful approach to the mapping of large areas in a semi-automated manner.In this paper, we introduce the study area, then outline the datasets and methods used (for land cover, palaeo-hydrology, and identifying anthropogenic signatures), and then provide some local examples illustrating these human-landscape relationships, based on detailed numerical analysis and modeling from remote sensing data coupled with contextual archaeological data.

Study Area
The study area includes the southeasternmost sector of the plateau of the Tassili n'Ajjer (Algeria), the valley of the Wadi Tanezzuft (Algeria/Libya), and a small part of the Tadrart Acacus (Libya), covering c. 3000 km 2 in total (Figure 2; see also par.2.2).Within this region, the Tassili is a rugged and largely inaccessible plateau [42], while the Tadrart Acacus features flat-topped structural terraces [10].Both areas are incised by wadis that have been the foci of past and current human occupations.The valley of the Tanezzuft, a 200 km-long strip running N-S and marking the international border between Algeria and Libya, separates the SE part of the Tassili from the SW sector of the Tadrart Acacus.The Tassili and the Tadrart were inscribed as UNESCO World Heritage Sites for their rock art in 1982 and 1985, respectively [43,44].The archaeology of the southeastern Tassili n'Ajjer plateau is still poorly known, and its restricted military zone status due to its proximity to the international border has hampered archaeological research in recent decades.However, important archaeological remains have been previously reported [43,45] and are clearly visible on satellite imagery [46].

Study Area
The study area includes the southeasternmost sector of the plateau of the Tassili n'Ajjer (Algeria), the valley of the Wadi Tanezzuft (Algeria/Libya), and a small part of the Tadrart Acacus (Libya), covering c. 3000 km 2 in total (Figure 2; see also par.2.2).Within this region, the Tassili is a rugged and largely inaccessible plateau [42], while the Tadrart Acacus features flat-topped structural terraces [10].Both areas are incised by wadis that have been the foci of past and current human occupations.The valley of the Tanezzuft, a 200 km-long strip running N-S and marking the international border between Algeria and Libya, separates the SE part of the Tassili from the SW sector of the Tadrart Acacus.The Tassili and the Tadrart were inscribed as UNESCO World Heritage Sites for their rock art in 1982 and 1985, respectively [43,44].The archaeology of the southeastern Tassili n'Ajjer plateau is still poorly known, and its restricted military zone status due to its proximity to the international border has hampered archaeological research in recent decades.However, important archaeological remains have been previously reported [43,45] and are clearly visible on satellite imagery [46].

SRTM and ASTER DEM
The SRTM 1 arc-second DEM topographic data (v.3) provides a horizontal resolution of 30 m and vertical accuracy of ±14 m.NASA released data covering regions outside the USA at this resolution, starting from September 2014.This near-global DEM was created using interferometric data acquired through the SRTM flown aboard the space shuttle Endeavor in February 2000 [47].A total of 16 SRTM tiles was used, covering a region slightly larger than the study area, to allow for the correct calculation of drainages prior to cropping to the exact study area, and were then mosaicked prior to processing.
The Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model Version 2 (GDEM V2) was released in October 2011 by the joint mission of the Ministry of Economy, Trade, and Industry (METI) of Japan and NASA [48].This product was generated using stereo-pair images collected by the ASTER instrument onboard Terra.The GDEM is provided at a 1 arc-second (30 m) resolution.The vertical accuracy of the DEM data is 17 m with 95% confidence without ground control point (GCP) correction for individual scenes.
Whilst in general SRTM data are considered to be more spatially consistent than ASTER since they are not affected by clouds, both systems show greater spatial and elevational error in very steep and rugged terrains, with ASTER having some advantage due to its nadir view.Both sensors have difficulties in smooth, flat terrain such as desert sand sheets, where little of the SRTM radar signal is reflected back to space, and where ASTER stereoscopy is hindered by the lack of earth's surface patterns to correlate the stereo views [49].Given these limitations, which may impact on the ability of either dataset to extract digital hydrological networks in desert regions, both SRTM and ASTER data were used in this project, and were compared for their overall performance.

Worldview 1-3
WorldView (WV) 1 to 3 data, acquired through grant aid from the DigitalGlobe Foundation (www.digitalglobefoundation.org), were used for mapping of archaeological features in the region.The area covered by WV2 imagery was used as a pilot for the creation of a vegetation map of the whole study area.WV2 covers, in its NE sector, a corner of Libyan territory where previous archaeological and geomorphological research has been carried out for the last 50 years, and also by one of the authors (e.g., 9).Therefore, we are confident that remotely sensed data from this region can be corroborated by field data.The whole dataset consisted of 3081.75 km 2 (Figure 2) of WV1-3 tiles.Launched between 2007 and 2014, the WorldView satellites have high spatial and increasingly higher spectral resolution sensors, featuring 0.50 m (WV1), 0.46 m (WV2) and 0.31 m (WV3) panchromatic imagery, 2 m multispectral imagery in 8 bands (WV2), 1.24 m multispectral imagery in 8 bands, and 8 short wave infrared bands at 3.7 m spatial resolution (WV3).The images were acquired in different dates, as reported in Table 1, under clear sky conditions.The acquired images, provided by DigitalGlobe™ at a Standard 2A level, had already been orthorectified and projected (Universal Transverse Mercator: UTM zone 32 North and WGS-84 Geodetic datum) and were not further processed geometrically by the authors.The WV2 images used for the Random Forest (RF, see Section 2.4.2) classification were radiometrically calibrated to obtain reflectance from radiance values, and were atmospherically corrected using the FLAASH (Fast Line-of-sight Atmospheric Analysis of Hypercubes) procedure in ENVI (Environment for Visualizing Images) v.5.3 software [50].FLAASH is a first-principles atmospheric correction tool that corrects wavelengths in the visible through near-infrared and shortwave infrared regions up to 3 µm, and incorporates the MODTRAN ® (MODerate resolution atmospheric TRANsmission) radiation transfer code to model atmospheres and aerosol types [51].It also corrects for the effects of atmospheric propagation on measurements acquired by air and space-borne systems.These corrections allow for the retrieval of accurate reflectance spectra.The atmospheric correction was performed using standard parameters with an average ground elevation of 1300 m, a Tropical atmospheric model, and no aerosol retrieval.

Delineation of Surface and Near-Surface Drainage Networks
Standard GIS hydrology procedures based on the D8 approach developed by Jenson and Domingue [52] were applied in ArcGIS v.10.3 to the mosaicked and void-filled 1 arc-second SRTM and ASTER data.As a first step, a fill procedure was applied to the DEMs in order to locate and fill spurious sinks and to ensure proper delineation of basins and streams.In fact, although sinks may occur naturally in a landscape (e.g., lakes, floodplains, ponds), in a derived DEM they are often manifested as errors due to the resolution of the data or by rounding of elevations to the nearest integer value, and do not necessarily correspond with actual features in the terrain [53].The sink fill function in ArcGIS progressively fills all sinks in the DEM by increasing their elevation value, and it iterates until all sinks within the specified z limit are filled.The z limit for sink filling was not specified, but was iterated until all sinks were filled, independent of the fill depth.Flow direction was determined by finding the direction of steepest descent from each cell with respect to surrounding cells.Each cell in the flow direction layer has its flow direction coded according to the D8 model in ArcGIS.Finally, to derive the stream network, a flow accumulation function was applied which calculates the number of upslope cells flowing to a particular location.Different thresholds (of 10,000 to 250 cells) were used to extract different raster stream networks from the flow accumulation file, and to allow examination of the hydrological characteristics of the region at a range of scales.In this manner, both small streams and major drainage systems can be differentiated.
After the raster stream network was vectorised, the results were verified through visual inspection against WV imagery where possible, and using high-resolution colour imagery in GoogleEarth (for the areas of the network that extend beyond the available WV imagery).This integrated approach allowed for visual interpretation of the results and the removal of spurious channels through pruning, in particular around paleolakes and natural depressions.

Defining Land Cover Classes and Reference Data Collection
Ground reference data could not be collected for this work due to the inaccessibility of the study area.Nevertheless the limited number of land cover classes that characterize the Sahara, and the high resolution of the multispectral and panchromatic imagery available for the study, allowed for the use of WV imagery itself for the choice of the training and testing classes.These broad classes were aimed at providing a general description of land cover, with particular focus on identifying different types of vegetation to contrast to other land cover types, such as sand and rock outcrops.A class identifying shadows in the rocky landscape (present in the proximity of steep slopes) was also created, since initial optimization tests on the dataset had shown that shadows played a significant role in confusing other classes.Seven broad classes were identified, and reference points were generated using classic heads-up techniques incorporating the colour (spectral), shape, texture, shading and pattern information of the available WV2 data, in a pansharpened image that combined panchromatic and multispectral information of the area examined.The reference points were then overlain on the WV2 image to create regions of interest (ROIs) to train and validate the classifier (Table 2) by randomly splitting the reference data into 70% training and 30% independent validation datasets (Table 2), which is a standard approach used for accuracy assessment in remote sensing studies [54].Decision learning trees such as classification and regression trees (CART) are commonly used for data mining and have been one of the most successful methods for undertaking supervised classification [55].To improve the accuracy of CART, Breiman [56] developed an ensemble learning technique called Random Forest (RF) by introducing the idea of bagging (bootstrap aggregating) to the decision trees.This method involves combining multiple decision trees where each tree contributes a single vote for the assignment of the most frequent class to the input data.Many binary classification trees (ntree) are built by RF using several bootstrap samples with replacements drawn from the original observations.Samples not in this bootstrap sample are called out-of-bag (OOB) samples.These OOB samples, which are about a third of the total dataset available, can be used to estimate the misclassification error and to measure the importance of each variable in the final model [56,57].A given number of input variables (mtry) at each node was randomly chosen from a random subset of the features, and the best split was calculated by utilizing only this subset of features.No pruning was performed and all trees in the forest are maximally-grown trees so as to ensure low bias [58].Mtry in this study is defined as the square root of the total number of spectral bands.In order to improve the classification accuracy, RF parameters (i.e., mtry and ntree) have to be optimized [56,59] and the default number of trees (ntree) is 500, while the default value for the number of variables (mtry) is the square root of the total number of spectral bands used in the study [56].A 10-fold grid-search approach based on the OOB estimate of error was used in this study to find the optimal combination for these two parameters, with the mtry value being varied from 1 to 5, and the ntree parameter varied from 500 to 10,000.The Image RF tool in EnMAP-Box was used to perform the RF classification.
In the present study, the package [60] developed for R statistical computing language was use to perform both the optimization and the classification algorithms.

Accuracy Assessment
The accuracy of the RF classifier was assessed using the independent test dataset (30%).A confusion matrix was constructed to compute the overall accuracy (OA), user's accuracy (UA) and producer's accuracy (PA) as criteria for evaluating the ability of the RF classifier to detect the classes present in the analysed landscape.OA is a ratio (%) of the number of correctly-classified samples against the number of test samples, while UA represents the likelihood that a sample belongs to a specific class and that the classifier accurately assigns it this class.PA expresses the probability of a certain class being correctly recognized.In addition, the kappa coefficient was calculated to provide a measure of the difference between the actual agreement, the reference data, and the classifier used to perform the classification, versus the likelihood of agreement between the reference data and a random classifier [54].

Anthropogenic Signatures in the Landscape
All of the panchromatic tiles of WV1-3 imagery were loaded in ArcGIS v.10.3.The study area was overlain with a grid consisting of 300 west-east strips, each of 250 m width.Each strip underwent systematic visual inspection, and every identified archaeological feature was digitised as a point located in the middle of the feature, and assigned to a thematic layer.All records were classified according to their morphological characteristics and size.Anthropogenic features that can be identified through imagery visual inspection mainly consist of tombs, monuments, settlements (including seasonal tented encampments and huts), forts (qsur), irrigation channels (foggara), crop fields, and stone enclosures of uncertain use (see Table 3).Other types of sites with an archaeological record, such as caves and rock shelters, cannot be identified through remote sensing data but they can be spatially located if they have been previously reported in the literature.Three issues are important in remotely sensed archaeology: (i) the intrinsically-problematic notion of defining an "archaeological site"; (ii) the correct identification of the archaeological evidence; and (iii) its correct chronological attribution.Archaeological landscapes are horizontal and vertical palimpsests, often lacking physical discontinuities.The adoption of the notion of "site" (issue i) has proven problematic in North African archaeology (for recent comments, see [1,61]) and in general [62,63].There are epistemological, methodological, and practical issues around the identification of "sites", especially where landscapes show a wide variety of different types of anthropogenic signatures.For this reason, we adopt the more neutral term "anthropogenic feature".Most of the features identified consist of well-delimited single structures (e.g., tomb, hilltop fort), although some more enigmatic (e.g., enclosures) or complex or composite features (e.g., modern pastoral campsite with a number of features within it) are less easily defined.To minimize the risk of misidentification (issue ii), we rely on operators' knowledge acquired through fieldwork in neighbouring areas of the SW Fazzan [1,9,29], Wadi ash-Shati region [5,64,65], and Fadnoun range [66].Direct experience of comparable archaeological evidence, coupled with information from published data, reduces the risk of misidentification and incorrect chronological attribution (issue iii), especially in view of the relatively high standardization of the majority of Saharan monuments and features.Twenty minutes of uninterrupted inspection have proven to be an average time for completing a single strip.The digitization of the entire study area took ca. 100 h.Verification and validation of the data were performed by randomly repeating the digitization of the same strip by different operators, in order to test for misidentification and/or false positives.Final identification of evidence, in particular of traces classified as "uncertain", is still undergoing revisions and subject to minor changes.

Statistical Analysis
In order to evaluate whether variation in the spatial density of anthropogenic features was determined at least in part by vegetation and drainage, a series of Poisson point process models have been fitted to the data and compared using the Akaike Information Criteria (AIC) and McFadden's pseudo-R 2 using the spatstat package [67,68] in R statistical language [60].Models were based on three covariates: (1) an elevation map based on the ASTER DEM; (2) an Euclidean distance map with respect to the drainage network; and (3) a 1-km bandwidth Gaussian kernel density map of the vegetation, as extracted using the Random Forest classifier.

Hydrological Network
We have produced a drainage network that optimizes the results of extracted ASTER and SRTM data at a range of scales.This network increases the accuracy of the available datasets (HydroSHEDS) for the study area (Table 3).The generated data successfully delineate the major and minor drainages and paleodrainage courses (Figure 3).This precision allows for high-detail regional and site-specific paleohydrological research, in combination with archaeological data, and other proxies identified in this study.
respect to the drainage network; and (3) a 1-km bandwidth Gaussian kernel density map of the vegetation, as extracted using the Random Forest classifier.

Hydrological Network
We have produced a drainage network that optimizes the results of extracted ASTER and SRTM data at a range of scales.This network increases the accuracy of the available datasets (HydroSHEDS) for the study area (Table 3).The generated data successfully delineate the major and minor drainages and paleodrainage courses (Figure 3).This precision allows for high-detail regional and site-specific paleohydrological research, in combination with archaeological data, and other proxies identified in this study.

Random Forest Parameters Setting
Random forest input parameters (ntree and mtry) were optimized using a grid search and a 10-fold cross validation (CV) method.The results illustrate that ntree values of 1000, 2500, 5500, 6500 and 7500 combined with mtry values of 2 and 3 produced the lowest OOB error (12.5%).The highest OOB error was produced by the highest mtry values of 8 (Figure 4).

Random Forest Parameters Setting
Random forest input parameters (ntree and mtry) were optimized using a grid search and a 10-fold cross validation (CV) method.The results illustrate that ntree values of 1000, 2500, 5500, 6500 and 7500 combined with mtry values of 2 and 3 produced the lowest OOB error (12.5%).The highest OOB error was produced by the highest mtry values of 8 (Figure 4).

Vegetation and Other LULC Classification
A land use/land cover (LULC) map was generated using a RF classifier to discriminate between cultivated fields, wooded vegetation (mainly Acacia trees) and shrubs from other land cover and land use types (Figure 5).
The RF classifier also provides a measure of the contribution of every variable (WV2 bands) in the LULC mapping as part of the classification process (Figure 6) through the mean decrease in accuracy measure.The mean decrease in accuracy is determined during the OOB error calculation.The more the accuracy of the RF decreases due to the exclusion (or permutation) of a single variable, the more important that variable is deemed, and therefore variables with a large mean decrease in accuracy are more important for classification of the data.This measure, illustrated in Figure 5, clearly indicates that the near-infrared bands (bands 7 and 8) were the most influential in the classification process since they present the highest mean decrease in accuracy values.The overall accuracy of the LULC classification drops by about 4.5% when the near-infrared bands are omitted from the classification model.

Vegetation and Other LULC Classification
A land use/land cover (LULC) map was generated using a RF classifier to discriminate between cultivated fields, wooded vegetation (mainly Acacia trees) and shrubs from other land cover and land use types (Figure 5).
The RF classifier also provides a measure of the contribution of every variable (WV2 bands) in the LULC mapping as part of the classification process (Figure 6) through the mean decrease in accuracy measure.The mean decrease in accuracy is determined during the OOB error calculation.The more the accuracy of the RF decreases due to the exclusion (or permutation) of a single variable, the more important that variable is deemed, and therefore variables with a large mean decrease in accuracy are more important for classification of the data.This measure, illustrated in Figure 5, clearly indicates that the near-infrared bands (bands 7 and 8) were the most influential in the classification process since they present the highest mean decrease in accuracy values.The overall accuracy of the LULC classification drops by about 4.5% when the near-infrared bands are omitted from the classification model.

Accuracy Assessment
Accuracy assessment was carried out in order to evaluate the prediction performance of the RF classifier using the independent test dataset.Table 4 describes the confusion matrix of the LULC types.The classification model has an overall accuracy of 82.01% with a Kappa value of 0.789.Table 4. Confusion matrix and classification of the accuracies of the seven classes of wooded vegetation (WV), crops (CR), rocks (RK), sand dunes (SD), shrubs (SH), shadow (SW), and water bodies (WB).The accuracies include the overall classification accuracy (OAC), kappa statistic, producer accuracy (PA) and user accuracy (UA).In general, all LULC classes achieved over 76% user's accuracy, with the exception of shrubs (SH), which had 70.27% due to spectral confusion with wooded vegetation (mainly Acacia trees), crops, rocks and sand dunes.All LULC classes achieved over 80% producer's accuracy, with the exception of crops (76.67%) and shrubs (59.09%).Crops were confused with wooded vegetation and shrubs, possibly due to a similar spectral response (according to the classifier).Shrubs were confused with wooded vegetation, crops and sand, which resulted in an overestimation of the

Accuracy Assessment
Accuracy assessment was carried out in order to evaluate the prediction performance of the RF classifier using the independent test dataset.Table 4 describes the confusion matrix of the LULC types.The classification model has an overall accuracy of 82.01% with a Kappa value of 0.789.Table 4. Confusion matrix and classification of the accuracies of the seven classes of wooded vegetation (WV), crops (CR), rocks (RK), sand dunes (SD), shrubs (SH), shadow (SW), and water bodies (WB).The accuracies include the overall classification accuracy (OAC), kappa statistic, producer accuracy (PA) and user accuracy (UA).In general, all LULC classes achieved over 76% user's accuracy, with the exception of shrubs (SH), which had 70.27% due to spectral confusion with wooded vegetation (mainly Acacia trees), crops, rocks and sand dunes.All LULC classes achieved over 80% producer's accuracy, with the exception of crops (76.67%) and shrubs (59.09%).Crops were confused with wooded vegetation and shrubs, possibly due to a similar spectral response (according to the classifier).Shrubs were confused with wooded vegetation, crops and sand, which resulted in an overestimation of the presence of vegetation in sand dune areas (easily confirmed by visual inspection of the classified image with WV2 panchromatic imagery covering the same area).This is most probably due to the morphology of the shrubs' canopy, which is not dense and permits high visibility of the sandy surface below.The classifier could therefore have confused these mixed spectral signatures (Figure 7).This confusion between sand dunes and shrubs not only reduces both user's and producer's accuracies, but contributes negatively to the overall combined accuracy.In arid areas the reflection of soils, rocky outcrops, and sand therefore poses serious challenges to the discrimination of sparse and patchy vegetation in these landscapes [69].

WV
presence of vegetation in sand dune areas (easily confirmed by visual inspection of the classified image with WV2 panchromatic imagery covering the same area).This is most probably due to the morphology of the shrubs' canopy, which is not dense and permits high visibility of the sandy surface below.The classifier could therefore have confused these mixed spectral signatures (Figure 7).This confusion between sand dunes and shrubs not only reduces both user's and producer's accuracies, but contributes negatively to the overall combined accuracy.In arid areas the reflection of soils, rocky outcrops, and sand therefore poses serious challenges to the discrimination of sparse and patchy vegetation in these landscapes [69].

Anthropogenic Features
Visual inspection has identified more than 7000 anthropogenic features in the study area (Table 5, Figure 8).Almost all grave types known to occur in the Central Sahara have been identified in the study area.The majority of these graves (n = 5261 tombs) are round stone structures 3-5 m in diameter, or drum shaped or bazinas (Figure 9), which are likely of Garamantian age.These tombs are generally clustered in necropoleis, representing the most tangible evidence of historical human occupation in the Fazzan.The study area is close to the Garamantian necropolis of Fewet [70], and includes the published graves of Ti-n-Alkoum [45].In contrast to the Garamantian burials, conical tumuli of the Late Pastoral age are more isolated, as confirmed by recent fieldwork in the southern reaches of the Tadrart Acacus [71], and can be of large sizes.The west side of the Wadi Tanezuft appears to be a key area for likely Garamantian occupation, not far from present-day oases.On the other hand, other complex graves are widely dispersed, but some parts of the study area seem to be devoid of such evidence, especially throughout the Tadrart Acacus.
Six fortified/compound settlements have also been recorded thus far.These are characterized by a stone enceinte that delimits the perimeter of the settled area.Some of these forts are located on hilltops, consistent with those already noted in the Wadi el Ajal [21] and Wadi ash-Shati [5].These forts reportedly have a long history of occupation, commonly from historic to medieval times but occasionally to the modern period.More than 600 features have been labelled as "uncertain", roughly about 10% of the whole dataset.These are features that cannot be definitively classified and should be subjected to ground-truthing.

Anthropogenic Features
Visual inspection has identified more than 7000 anthropogenic features in the study area (Table 5, Figure 8).Almost all grave types known to occur in the Central Sahara have been identified in the study area.The majority of these graves (n = 5261 tombs) are round stone structures 3-5 m in diameter, or drum shaped or bazinas (Figure 9), which are likely of Garamantian age.These tombs are generally clustered in necropoleis, representing the most tangible evidence of historical human occupation in the Fazzan.The study area is close to the Garamantian necropolis of Fewet [70], and includes the published graves of Ti-n-Alkoum [45].In contrast to the Garamantian burials, conical tumuli of the Late Pastoral age are more isolated, as confirmed by recent fieldwork in the southern reaches of the Tadrart Acacus [71], and can be of large sizes.The west side of the Wadi Tanezuft appears to be a key area for likely Garamantian occupation, not far from present-day oases.On the other hand, other complex graves are widely dispersed, but some parts of the study area seem to be devoid of such evidence, especially throughout the Tadrart Acacus.
Six fortified/compound settlements have also been recorded thus far.These are characterized by a stone enceinte that delimits the perimeter of the settled area.Some of these forts are located on hilltops, consistent with those already noted in the Wadi el Ajal [21] and Wadi ash-Shati [5].These forts reportedly have a long history of occupation, commonly from historic to medieval times but occasionally to the modern period.More than 600 features have been labelled as "uncertain", roughly about 10% of the whole dataset.These are features that cannot be definitively classified and should be subjected to ground-truthing.

Point Process Models
When environmental factors are examined individually, elevation showed the best fit to data, with a McFadden's pseudo-R 2 of 0.1652 (Table 6).Drainage and vegetation both performed poorly individually, but models combining both variables with elevation (either considering only the main effects or including interaction) showed a drastic decrease in the AIC, and a pseudo-R 2 over 0.2, suggesting a fairly good fit to the data (see [70]: 306-307 for interpretation of pseudo-R 2 in relation to coefficients of determination).

Point Process Models
When environmental factors are examined individually, elevation showed the best fit to data, with a McFadden's pseudo-R 2 of 0.1652 (Table 6).Drainage and vegetation both performed poorly individually, but models combining both variables with elevation (either considering only the main effects or including interaction) showed a drastic decrease in the AIC, and a pseudo-R 2 over 0.2, suggesting a fairly good fit to the data (see [70]: 306-307 for interpretation of pseudo-R 2 in relation to coefficients of determination).

Discussion
High-resolution imagery, such as the WorldView data we have used in this study, is undoubtedly a powerful tool for identifying, interpreting, and analysing the late Holocene record of the central Sahara.The high number of anthropogenic features recorded, along with their identification, clearly demonstrates the significance of Earth Observation techniques for the mapping and monitoring of cultural heritage.Previous direct knowledge of the area is, however, of added value to the remote sensing data in terms of field validation and in providing radiometric dating control.In the absence of previous fieldwork in the study area, it is likely that the identification and interpretation of anthropogenic features would have been more difficult.
In addition, our integrated research approach is a useful way to explore and interpret cultural heritage in regional environmental context(s).Our analyses confirm the relevance of vegetation and drainage for understanding spatial variations in archaeological evidence.The fitted point process model suggests that the interplay of vegetation, drainage and elevation successfully explains the global variability in the distribution of anthropogenic features (Table 6).However, different types of vegetation, drainage or different ages of anthropogenic features were not disaggregated for this study, nor different types of local and inherent spatial dependence (e.g., clustering or repulsion) were not considered in the point process model [72][73][74].Furthermore, the model so far does not include the limited number of rock art sites published (e.g., [11]) covering the last 5000 years.
Although full results from the wider study area are incomplete, the examples presented here illustrate the predictive and interpretive potential of combining present-day vegetation and hydrological data.This study also confirms that classifying vegetation in desert environments is not trivial [69].Traditional pixel-based image methods, such as the ones used in this study, are in fact limited since image pixels are not true geographical objects and pixel topology is limited and overlooks the spatial photo-interpretive elements such as texture, context and shape [75,76].Moreover, with high spatial resolution imagery such as WV2, the variability implicit within the class sometime confuses the pixel-based classifiers and results in lower accuracy, in particular for vegetation classes.Therefore, we believe that misclassifications of the vegetation classes seen in this study can be attributed to the high spectral variation within the same vegetation class.To improve our results, alternative approaches that integrate object-based classification and advanced classifier algorithms should be used.Another possible way to improve our results could include the use of Linear Discriminant Analysis (LDA) [41], which can further analyze spectral signatures and thus may redefine vegetation classes not at pixel level.In the case of shrubs, whose spectral signature is a combination of vegetation and bare ground, segmentation processing (i.e., drawing polygons including a number of pixels) and analysis of segments could also be used instead of individual pixels [77].

Conclusions
This study shows that the use of satellite data at varying resolutions in arid lands can build a baseline for the application of spatial based analyses in a GIS environment.The research outlined here will provide further details on patterns of settlements and land use, and human-environment interactions in the late Holocene Sahara.The identification of (some of) the potential factors affecting the distribution of anthropogenic features in the area is a key result.The fitted point process model presented here showcases a small glimpse of the potential benefits of adopting similar (and more refined) spatial analyses for future work.There are two wider implications of our study: (1) an integrated research framework for remote sensing studies of desert landscapes has the potential to enhance our understanding of human-environment relations over millennial timescales.This methodology can be applied to other regions with long records of human occupation such as the Simpson Desert (Australia), Gobi Desert (China/Mongolia), Arizona/New Mexico (southern USA), Kalahari (southern Africa) and Jordan/Syria (Middle East); (2) The use of remote sensing allows for cost-effective research to continue in areas that are geopolitically unstable, and the approach can be used as an effective monitoring tool for sites of archaeological or built heritage value.The role of remote sensing in this context has been forcibly shown by images of the recent destruction of World Heritage Sites in Syria, such as Palmyra.The use of remote sensing and environmental modeling approaches may be able to capture data from other sites similarly under threat.

Figure 1 .
Figure 1.The schematic GIS workflow used in this study.Blue boxes include types of satellite imagery.Processing steps are in ovals, representing ArcGIS processing (yellow), ENVI processing (orange) and R processing (purple).Red boxes feature the output datasets.The green box represents the final output.

Figure 1 .
Figure 1.The schematic GIS workflow used in this study.Blue boxes include types of satellite imagery.Processing steps are in ovals, representing ArcGIS processing (yellow), ENVI processing (orange) and R processing (purple).Red boxes feature the output datasets.The green box represents the final output.

Figure 2 .
Figure 2. Details of the study area.(a) Geopolitical setting (base map: ESRI); (b) Digital Elevation Model (DEM) at a spatial resolution of 1 arc-second (30 m) overlain with the area of interest (black polygon, which includes parts of the Tassili n'Ajjer and the Tadrart Acacus).SRTM data courtesy of the United States Geological Survey (USGS, 2014, base map: SRTM v4.1, www.cgiar-csi.org);(c) detail of the study area boxed in (b) featuring the WorldView (WV) tiles used for mapping purposes (blue: WV1, red: WV2, green: WV3); for the characteristics of the WV imagery, see Section 2.2.2.

Figure 2 .
Figure 2. Details of the study area.(a) Geopolitical setting (base map: ESRI); (b) Digital Elevation Model (DEM) at a spatial resolution of 1 arc-second (30 m) overlain with the area of interest (black polygon, which includes parts of the Tassili n'Ajjer and the Tadrart Acacus).SRTM data courtesy of the United States Geological Survey (USGS, 2014, base map: SRTM v4.1, www.cgiar-csi.org);(c) detail of the study area boxed in (b) featuring the WorldView (WV) tiles used for mapping purposes (blue: WV1, red: WV2, green: WV3); for the characteristics of the WV imagery, see Section 2.2.2.

Figure 3 .
Figure 3.Comparison of the HydroSHEDS hydrological product (yellow) with the hydrological network derived from the ASTER DEM (using 1000 cells) (blue) used in this study.

Figure 3 .
Figure 3.Comparison of the HydroSHEDS hydrological product (yellow) with the hydrological network derived from the ASTER DEM (using 1000 cells) (blue) used in this study.

Figure 4 .
Figure 4. Results of random forest optimization.The OOB error was calculated using a 10-fold cross validation (CV) and the training data sets.

Figure 4 .
Figure 4. Results of random forest optimization.The OOB error was calculated using a 10-fold cross validation (CV) and the training data sets.

Figure 5 .
Figure 5. Land use/land cover (LULC) map generated using a RF classifier (Satellite image courtesy of the DigitalGlobe Foundation).

Figure 5 .
Figure 5. Land use/land cover (LULC) map generated using a RF classifier (Satellite image courtesy of the DigitalGlobe Foundation).

Figure 6 .
Figure 6.Illustration of the importance of WV2 bands in vegetation and other LULC classifications, calculated using the mean decrease in accuracy index.The highest mean decrease in accuracy indicates the most important bands.Bands 7 and 8 are therefore considered to be most important and Band 5 least important.

Figure 6 .
Figure 6.Illustration of the importance of WV2 bands in vegetation and other LULC classifications, calculated using the mean decrease in accuracy index.The highest mean decrease in accuracy indicates the most important bands.Bands 7 and 8 are therefore considered to be most important and Band 5 least important.

Figure 7 .
Figure 7. Spectral signatures of the RF classifier, before (a) and after (b) atmospheric corrections.

Figure 7 .
Figure 7. Spectral signatures of the RF classifier, before (a) and after (b) atmospheric corrections.

Figure 8 .
Figure 8. Map of anthropogenic features recorded in the study area (boxed) (satellite image courtesy of the DigitalGlobe Foundation).

Figure 8 .
Figure 8. Map of anthropogenic features recorded in the study area (boxed) (satellite image courtesy of the DigitalGlobe Foundation).

Figure 9 .
Figure 9. Selected examples of anthropogenic features recorded in the study area.(a) Garamantian necropolis and hilltop qasr; (b-d) different stone cairns from the Late/Final Pastoral period; (e) enigmatic stone enclosures; (f) Tuareg campsite (panchromatic satellite imagery courtesy of the DigitalGlobe Foundation).

Figure 9 .
Figure 9. Selected examples of anthropogenic features recorded in the study area.(a) Garamantian necropolis and hilltop qasr; (b-d) different stone cairns from the Late/Final Pastoral period; (e) enigmatic stone enclosures; (f) Tuareg campsite (panchromatic satellite imagery courtesy of the DigitalGlobe Foundation).

Table 1 .
Synoptic characteristics of WorldView imagery used in this study.

Table 2 .
Number of training and validation datasets collected for different land cover classes in the study area.

Table 3 .
Comparison of the extracted hydrological data using different datasets in this study (ASTER, SRTM), with existing data (HydroSHEDS).

Table 5 .
Summary of different anthropogenic features identified in this study.

Table 5 .
Summary of different anthropogenic features identified in this study.Approx.Chronology N Complex Graves Late/Final Pastoral (c.3000-1000BC) Antenna 9 Crater 13 Crater on platform 2