Potentialities and Limitations of Research on VHRS Data: Alexander the Great’s Military Camp at Gaugamela on the Navkur Plain in Kurdish Iraq as a Test Case

Abstract

This paper presents a selected aspect of research conducted within the Gaugamela Project, which seeks to finally identify the location of one of the most important ancient battles: the Battle of Gaugamela (331 BCE). The aim of this study was to discover material remains of the Macedonian military camp on the Navkur Plain in Kurdish Iraq. For this purpose, three very high resolution satellite (VHRS) datasets from Pleiades and WorldView-2 were acquired and subjected to multi-variant image processing (development of different color composites, integration of multispectral and panchromatic images, use of principle component analysis transformation, use of vegetation indices). Documentation of photointerpretation was carried out through the vectorization of features/areas. Due to the character of the sought-after artifacts (remnants of a large enclosure), features were categorized into two types: linear features and areal features. As a result, 19 linear features and 2 areal features were found in the study area of the Mahad hills. However, only a few features fulfilled the expected geometric criteria (layout and size) and were subjected to field groundtruthing, which ended in negative results. It is concluded that no traces have been found that could be interpreted as remnants of an earthen enclosure capable of accommodating around 47,000 soldiers. Further research perspectives are also suggested.

1. Introduction

1.1. Historical Framework

The Battle of Gaugamela, which occurred in 331 BCE between the Macedonian and Greek forces commanded by Alexander the Great and the Persian forces under Darius III, respectively, has often been referred to as one of the most significant battles in the history of the Mediterranean world [1,2,3,4,5,6,7]. Indeed, its final result (Alexander’s sweeping victory) led to both the effective collapse of the Persian Empire—an empire that spanned two centuries—and the emergence of a new age, now commonly labeled the Hellenistic period. It was the Hellenistic period that brought about the unprecedented exchange of Greek and Oriental cultures throughout the ancient Near East (as far as modern Pakistan) and as such laid the foundations for the cultural heritage of the Mediterranean and Middle Eastern worlds that continues to this day (see Figure 1A).

Despite the great importance of this battle, its exact location is disputed. This situation results from the state of the extant historical sources [8,9,10,11,12,13,14]. First, the extant ancient sources do not provide us with very precise geographical and topographical information (at least not to an extent that would satisfy modern geographers and cartographers), as they either take such information for granted or are simply not interested in such practical aspects. Second, some important pieces of information that can be gleaned from various ancient accounts actually contradict each other. For instance, according to Arrian (Arr. 3.15.5, 6.11.4–6), the battlefield was located around 600 or 500 (two slightly different traditions quoted by Arrian) stades from Arbela (modern Erbil), which equals ca. 111 or 92 km. In contrast, according to Curtius Rufus (Cur. 4.9.9–10), the battlefield was 60 stades away from the Lykos River (at the point where the Persians crossed it), which equals ca. 14 km. Given the fact that the three most frequently used modern crossings of the Lykos River (that is, the modern Great Zab River) are located some 40.8 km (Chammah), 38.9 (Eski Kelek), and 57.3 km (Al Kuwayr) away from Erbil, respectively (only to give a parallel!), the distances given by Arrian and Curtius exclude each other—Curtius locates Gaugamela much closer to Arbela than does Arrian.

As a result, it comes to no surprise that modern scholars have not agreed on one location for the Gaugamela battlefield. Thus, several locations have been suggested. In fact, all previous identifications go back to early European travelers and scholars from the eighteenth to the early twentieth centuries. Several sites have been suggested to date as the location of ancient Gaugamela: Karamleis [15], Qaraqosh [2], Tell Aswad [16], (a mound south of) Wardak [17], and Tell Gomel [18,19,20]. All of these identifications have recently been conveniently labeled as the southern (including Karamleis, Qaraqosh, Tell Aswad, and Wardak) and northern (Tell Gomel) hypotheses [21]. In fact, on the authority of A. Stein, Karamleis, in particular, has become the most well-known identification of Gaugamela, south of Jebel Maqlub [22,23]. In turn, Tell Gomel is the only site north of Jebel Maqlub that enjoys an association with the famous battlefield. As a result, scholars have chosen to side with only one of these identifications for many decades—either the Tell Gomel area [3,4,24,25,26,27,28,29,30,31,32,33,34,35,36] or the Karamleis region [5,37,38,39,40] (see Figure 1B).

In contrast, the recent period of stability in the Kurdish autonomous region in northern Iraq has brought modern archaeology into full bloom [41,42,43]. Many international teams commenced their work in Kurdish Iraq, starting in 2011, with new teams arriving every year. One of these teams is the Land of Nineveh Archaeological Project (LoNAP hereinafter), which began its work in 2012 [44,45,46]. LoNAP is a wide-ranging multidisciplinary research project carried out by the Italian Archaeological Mission in Assyria (IAMA), and its aim is to study the archaeological landscape of the region of Dohuk and record, conserve, and promote the rich cultural heritage of this part of Iraqi Kurdistan. The project investigates a 3000 sq. km area by means of archaeological survey and excavations (see Figure 1B). The license of LoNAP includes the site and vicinity of Gir-e Gomel (Tell Gomel), which is one of the two major identifications of Gaugamela. Against this background, the Gaugamela Project was created in 2018 (preceded by a dedicated survey by LoNAP in 2016), and its specific aim is to establish the exact location of the Battle of Gaugamela (331 BCE) using a well-known approach combining the methods of ancient history and landscape archaeology (use of ancient textual evidence, scholarly literature, past and recent cartographic data, satellite remote sensing imagery, geographic information systems capabilities, and fieldwork). As of 2021, the Gaugamela Project is ongoing; closure is expected around 2022.

One of many questions relevant to the identification of the exact location of the Battle of Gaugamela is the possible existence of the material remains of an ancient battlefield. This may be a very problematic issue, given that over 2000 years have passed since the battle [47,48]. In this context, the possibility has been considered that such traces may be buried under the current surface and could be revealed with the help of modern technologies, especially remote sensing. One of such traces could be the Macedonian military camp set up next to the battlefield. Two ancient sources (Arrian, Anabasis of Alexander, 3.9.2–3, and Quintus Curtius Rufus, Histories of Alexander the Great, 4.12.15–24) claim that the two armies could not see each other because of intervening hills at a distance of 60 Greek stades (ca. 11 km). They gained a full view of each other once the Macedonians reached the hills 30 stades away from the Persian final positions (ca. 5.5 km). On the eve of battle, the Macedonians, who counted 47,000 soldiers, set up their final camp on these hills directly overlooking the battlefield. In the course of our topographical survey, the hills located west of the modern town of Mahad were suggested as the place of the final Macedonian camp (see Figure 2). This identification is based on several premises. First, the hills are located in the vicinity of Tell Gomel, the name of which goes back to the ancient name of Gaugamela [45]. Second, the hills are located along the assumed line of approach taken by Alexander’s troops from the modern Mosul Dam area to Tell Gomel [49]. The existence of a transportation route in this area is clearly conditioned by the local terrain (the Navkur Plain is sharply delimited by the line of the Zagros Mountains to the north and a range of hills to the south, including the prominent outcrop of Jebel Maqlub) and can also be concluded on the basis of the archaeological data recently gathered by LoNAP (major patterns of settlement and several remains of hollow ways in this area). Third, the location of the hills matches the details mentioned by the ancient sources with regard to distance and visibility: the hills are about 5–6 km distant from the Al-Khazir River (which must have marked the contour of the Persian line) and offer a commanding view of the entire plain toward Tell Gomel. At the same time, the hills form a total obstacle for visibility between any position located 5–6 km to the west (thus about 10–12 km west of Tell Gomel) and the assumed Persian position to the east (which has been verified both on the ground and using Geographic Information System (GIS) visibility tools).

How was the Macedonian camp constructed? The vocabulary employed in ancient sources for the description of the Macedonian camp indicates the use of a protective trench (taphros) and a palisade (charax) from stakes, but the palisade stakes must have been transported with the army and thus removed when breaking the camp [50,51,52]. However, no Macedonian marching camp has ever been archaeologically discovered. What remains to be done in such cases is to look for other analogies. It follows from ancient sources that the Macedonian camps were much different from Greek camps [53,54,55]. In contrast, it is sometimes acknowledged that the Macedonian marching camps were predecessors of the later Roman marching camps [52,53,56,57]. Indeed, both the Macedonian and Roman military marching camps followed the basic concept of a space enclosed by a trench and rampart with a palisade [58,59]. The trenches around the marching camp would not have been very deep or wide—about 1 or 2 m deep and wide—but together with the low wall behind the trench, there would have been a defensive barrier almost 3 m wide and 3 m high [60,61]. If the Macedonian enclosure was close to the Roman model, one would expect a regular shape—roughly square or rectangular [58,59]. However, a certain amount of adaptation of the regular shape to the existing features of the terrain was also practiced [59]. Roman analogies could also be considered when calculating the size of the marching camp. Nevertheless, different assessments of the camp size with regard to the number of troops can be found in the literature. For instance, a space of somewhere between 30 and 80 ha can be suggested (depending on the density inside the camp) for 47,000 soldiers if A. Richardson’s or S. Kayes’ research is to be followed [62,63,64]. In contrast, in light of Hanel’s data, 42.7 ha for 7000 cavalry and 144–192 ha for 40,000 foot soldiers would have been required [59]. At any rate, 47,000 Macedonian troops would indeed have required a considerable space for their camp.

The area in question has belonged to the license area of LoNAP since 2012 [44,45]. An archaeological survey in the LoNAP license area was preceded by the systematic examination of aerial and satellite imagery (especially CORONA images). This, in turn, led to the identification of many archaeological sites (about 1100), which were later verified on the ground by means of an archaeological survey. Given the huge size of the LoNAP study area (about 3000 km²), an extensive, mixed-survey strategy was implemented based on a motor vehicle survey combined with field walking. What is more, an intensive off-site field survey was also conducted in chosen locations, including the vicinity of Tell Gomel. Namely, the area of Tell Gomel was surveyed by means of transect walking in the 25 km² area around the site of Tell Gomel in 2015 and 2016. However, no direct material evidence related to the Battle of Gaugamela has been revealed in the course of LoNAP’s research. This fact convinced the authors that a new attempt must be based on new tools, and the question has been posed as to whether the remains of the Macedonian camp could be revealed with the use of very high-resolution, multispectral satellite images, as was the case with similar military embankments in Mesopotamia at Dura Europos and Hatra. Namely, the famous siegeworks at Hatra, dated to ca. 240/241 CE, were detected using a QuickBird image (0.6 m pixel) taken in October 2004 [65,66], while two previously unknown military enclosures at Dura Europos, dated to between 115 and 256 CE, were discovered using a QuickBird-2 image taken on 2 January 2006 [67].

1.2. Potentials and Problems of the Earth Observation Imagery

The large diversity of Earth observation imagery allows for the registration of information about terrain surface in various spatial resolutions and in a wide spectrum of electromagnetic radiation (visible, infrared, and microwave). This is accompanied by the increasing availability of image data, associated with the dynamically growing number of recording platforms and increasing resources for free archival data (for instance, USGS-Landsat and UE-Sentinel). At the same time, aerial and satellite photos, primarily in the visible range, are made available interactively on popular websites (Google Maps, Bing), and their use sometimes offers benefits for archaeology [68,69].

Despite the fact that there are many resources available for current satellite data, the progression of urbanization processes, growing acreage, and war damage make the time of the data’s registration the key issue in archaeological research in the Middle East [43,51] because these processes effectively mask or destroy historical artifacts [70]. In this respect, archival satellite images from the Cold War period (especially the CORONA program) are an extremely valuable source of data, and their value for archaeological research has already been proven a number of times [71,72,73]. In light of the fact that the use of the CORONA images by LoNAP has not revealed the suspected features, the decision has been made to employ current very high resolution satellite (VHRS) satellite data for this purpose. It has been assumed that the very high spatial resolution of the images, the use of spectral bands in the visible and infrared ranges, and the possibilities of digital image processing will allow for the effective use of data despite the problems mentioned above concerning the degradation of historical features that has been progressing in recent decades. This assumption is in line with the good results achieved with the use of the visible and near-infrared (NIR) range of the spectrum, especially the RedEdge band, in satellite archaeology [74].

Due to the problematic nature of the sought-after artifacts (imprecise historical descriptions and a lack of archaeological analogies), photointerpretation of image data with field groundtruthing was adopted as the basic research tool. Advanced methods of automated image data processing such as classifications of various types, including the use of neural networks or machine learning in particular, require a large knowledge base and reference data (a training set and a test set), which has not been available in the case of our research. The principles of photointerpretation, including the importance of direct (such as shape, color, and size) and indirect (such as associations and shade) photointerpretation features, as well as the use of color composites in false colors, have been widely discussed by authors such as Jensen [75] and Lillesand et al. [76]. A novelty in our research was the use of multi-variant imaging capabilities for the purpose of preparing additional materials for visual interpretation [76,77,78,79], including the use of decorrelation stretch, so far marginally used in archaeological research. In the subject literature, it is referred to as a procedure useful for enhancing and correcting a graphic reproduction of photographic information about individual archeological sites, such as ancient Egyptian rock paintings (ground registration) [80] or boundaries of non-extant buildings of an ancient Roman city (large-scale aerial photos) [81]. The novelty in our case has been that we tried to use it to detect—in a relatively large area—new, unclear archaeological features.

It is evident that materials and techniques employed in this study are only a small part of the possibilities offered by remote sensing [82] and employed in archaeological prospection [83,84]. However, due to the budget and time constraints of the project, other recording techniques such as radar, Ground Penetrating Radar (GPR), thermal imaging, hyperspectral imaging, and laser scanning were not employed. This type of data acquisition, assuming the need for their high spatial resolution, requires a low registration height (aerial, UAV) [81,85], resulting in large datasets for a relatively small imaging area. This translates into high registration costs and long processing times. Given these factors and the limited accessibility of the area due to the political situation, their use has been considered as an option in the project, but only at a later stage when the probability of the occurrence of sought-after features would be high. It is worth emphasizing that this type of data have proven to be highly useful in archaeological prospection. Synthetic aperture radar (SAR), using a satellite sensor in the microwave frequency range, can be an especially potent tool for archaeological spaceborne prospection [86]. For instance, SAR was successfully employed to detect the course of a Roman road from the second century BCE across wet, untilled agricultural fields in southern Spain [87]. It was also used to detect buried walls of the Roman fortress Qreiye at a depth of up to 25 cm in the dry desert soils of eastern Syria [88]. Likewise, it was also effectively used to detect buried archaeological features both in vegetated areas in the vicinity of Rome [89] and in desert areas in the North Sinai desert [90]. Both of these indicators—plant and soil—were crucial in the aforementioned study, so the use of this type of technique would be justified at the next, more detailed stage of our project. Similarly, the use of hyperspectral data has proven effective in archaeological prospection in the works of Cavalli et al. [91], Aqdus et al. [92], and Doneus et al. [93], among others. A promising option for the search of subsurface archaeological remains is the use of automatic classification of hyperspectral data. In the areas covered with vegetation, the most effective are the yellow-green and RedEdge bands; in areas with uncovered soils, the thermal bands are most effective [94]. A large number of different vegetation indices in the hyperspectral data have been examined by Agapiou, and it has been shown that next to the most widely used vegetation indices in satellite archaeology (especially the Normalized Difference Vegetation Index (NDVI), Spectral Response (SR) ), several other indices can assist in the detection of crop marks [94]. Likewise, the analysis of altitude models based on laser scanning is an effective tool mainly when searching for archaeological features in areas covered by other features, such as dense vegetation in tropical areas [95] or water [96]. Not surprisingly, the extensive use of Light Detection and Ranging (LiDAR) in one research project most recently led to the spectacular discovery of 66 Roman military camps in northwestern Spain [97]. The use of a high-quality Digital Elevation Model (DEM) in this study or the availability of stereograms would improve the initial verification of detected features, limiting subsequent field groundtruthing. Unfortunately, the DEM data available from global models for the area of analysis are limited in accuracy and their usefulness is questionable. Their use would require verification based on specially conducted additional research.

Free, open-access satellite data from medium-resolution missions such as Landsat TM/ETM/OLI or Sentinel-2 were not used either, as this kind of data were not especially useful for the project. These missions, which additionally contain mid-infrared bands, provide great interpretation capabilities at the regional and sub-regional levels, but they do not allow the detection of relatively small structures, anomalies, and narrow linear features. However, they can be a valuable source of information about general processes occurring in a wider window of time (even since the 1980s), affecting the general archaeological context (e.g., floods, fires, earthquakes, war damage, major fieldwork, urbanization, etc.). In this case, as good knowledge of the area was achieved by LoNAP over several field seasons, this type of analysis was unnecessary.

2. Materials and Methods

2.1. Data

The available archival databases were searched for the area of interest in June 2018 (https://www.intelligence-airbusds.com, http://iohs.euspaceimaging.com (accesssed date 11 June 2018)). Several criteria were defined: a maximum value of spatial resolution of pixels at 2 m for multispectral (MS) imagery, a maximum value of spatial resolution of pixels at 0.5 m for panchromatic (PAN) images, a maximum percentage of cloud cover at 10%, and an imaging angle below 20°. The choice of data was limited to the available very high resolution satellites (VHRSs) that had proven reliability and recorded multispectral images in a visible and a near-infrared range of the spectrum, e.g., Ikonos, QuickBird, GeoEye, Pleiades, WorldView, KompSat, and Spot 6/7. Eventually, three VHRS datasets, registered in 2013–2015, were chosen from Pleiades and WorldView-2 (see

Table 1. Characteristics of satellite data used in the analyses.

Satellite	Pleiades 1A, 1B © CNES (2013, 2015), Distribution Airbus DS	WorldView-2 DigitalGlobe
Repository	11 April 2013 DS_PHR1A_201304110806519 12.12+24.12.2015 DS_PHR1B_201512120807314 DS_PHR1B_201512240814454	28 September 2014 ID10300100388F8C00
Panchromatic (PAN) (micrometers)	0.47–0.83	0.45–0.80
Multispectral (MS) (micrometers)	0.43–0.55 0.50–0.62 0.59–0.71 0.74–0.94	0.40–0.45 (coastal) 0.45–0.51 0.51–0.58 0.58–0.62 (yellow) 0.63–0.69 0.70–0.75 (red edge) 0.77–0.90 near-infrared (NIR 1) 0.86–1.05 (NIR 2)
Resolution PAN/MS	0.5 m/2 m	0.5 m/2 m
Radiometric resolution	12 bits	11 bits

). Due to the lack of a reliable DEM for this area and the lack of detailed geodetic measurements, it was decided not to perform image orthorectification but rather to use already processed products—orthophotomaps with the WGS84-UTM38 reference system.

Both Pleiades and WorldView-2 record a PAN band with a spatial resolution four times higher than the multispectral image. Increasing the resolution of color composites from 2 m to 0.5 m significantly increases their interpretative values. In the case of Pleiades images, orthophotomaps with a resolution of 0.5 m were ordered; these were already integrated (MS + PAN) with an algorithm used by the company ASTRI (the algorithm is not disclosed). The synthetic image is characterized by high suitability for photointerpretation. Because the integrated product was ordered, the data were 4 bands with a resolution of 0.5 m (there were no source data, i.e., one PAN image with a resolution of 0.5 m and four MS images with a resolution of 2.0 m). The WorldView-2 data were obtained in the form of 8 spectral bands with a spatial resolution of 2m and 1 PAN band with a spatial resolution of 0.5 m.

The decision on the selection and diversity of data was made in light of other published studies and the experience of researchers from the area of interest. In the case of winter and spring dates, the dominant interpretative factors were differences in soil moisture and vegetation, proven as extremely useful indicators of archaeological features in previous studies on Mesopotamian landscapes [70,71,73]. The early autumn season was chosen as an additional contrasting material. The use of spectral bands in the present research, including those in the non-visible field, as well as the digital processing of this data, brings new interpretative possibilities, especially in the use of plant indicators [70]. The scope of the image data processing proposed here is clearly wider than simple operations performed on CORONA images in similar studies. Although procedures such as principle component analysis (PCA) transformation and image fusion are well known in remote sensing, the comprehensive use of all these elements for photointerpretation in the archaeological application of remote sensing is novel. Consequently, it has been hoped that their use may offer the chance for a synergistic effect, leading to the detection of archaeological traces buried under the current surface.

2.2. Data Processing for Image Enhancement: Methodology

In the case of the present study, automated methods of obtaining information from images were not taken into account due to the lack of unambiguous patterns of the features being sought after. Namely, the main target was the discovery of the location where the Macedonian military set up on the eve of the battle next to the battlefield. The problem, however, is that historical sources do not provide precise details about the form of the Macedonian camps (especially their layout and size). What is worse, no Macedonian camp has ever been discovered archaeologically, and consequently, our assessment of the layout and size of the Macedonian camp can only be hypothetical (and in fact is based on the later analogies of Roman marching camps). In this light, it has been decided to rely on the manual photointerpretation method, with computer support for image enhancement.

The methodology employed in this study enabled the search for artifacts, including structures and places where geometrized areas/elements (most often humidity and plant cover known widely as marks, especially soil and crop marks) present on the surface of the land could be indicative of events from over 2000 years ago and are not justified in relation to the current land use [97,98]. The complementary use of various image data may be of assistance in such cases, which should increase the possibility of detecting anomalies on the surface of the area of interest.

In areas where, for example, the geometrization of specific vegetation occurs, local changes in soil or water properties should be checked for potential traces of the sought-after features. In cases like these, indirect features of photointerpretation are used to associate the existence of vegetation, soil, or water conditions with traces of features that are no longer visible on the surface level. This task seems to be more difficult than even the search for hollow ways detected in the panchromatic data from the 1960s and 1970s, where direct features such as width, lighter tone (different from the environment), linear shape, and appropriate length allowed for their detection [72,99,100,101,102]. Given the aforementioned historical analogies of the Roman marching camps, it has been estimated in our research that the physical features of the Macedonian camp (earthen rampart about 1–2 m high and wide; trench adjacent and parallel to the rampart of the camp) would translate into the following characteristics of the remains: linear elements of unknown shape and course limiting the area of the camp, each side from 0.5 to 1.5 km long. Linear elements may also have been partly natural and not anthropogenic (because of the adaptation of the camp features to the existing features of the terrain). A detailed list of the photointerpretation features of the sought-after linear features with explanations is presented in

Table 2. Sought-after linear elements and related photointerpretation features anticipated in the satellite imagery.

Possible Current State	Photointerpretation Features, Associations, Indicators	Abbreviations in Figures
Partly leveled rampart and trench through erosion and/or human activities (agriculture). Geographic Information System (GIS) encoding linear vector Linear exist (LE)	Direct features: - exposed contour of the course of ramparts - greater reflection of a convex feature (rampart: better illuminated, more dry) than of a concave feature (trench: not as well illuminated, more humid) - texture change (cutting of existing textures on uncovered soil, dense vegetation) Indirect features, associations: - the shadow of the feature (rampart), the diversity of illumination of the rampart slopes, depending on the direction of sunlight - different vegetation associated with varied sunlight and humidity of the rampart and trench (as well as in relation to the environment) - possible accompaniment of geomorphological features, paths, streams	c (contour) s (spectral) t (texture) sh (shadow) v (vegetation) g (geomorphology) p (paths) w (watercourses)
Completely leveled rampart, silted up with natural processes and/or a trench covered by human activity (agriculture, urbanization) GIS encoding linear vector Linear non-exist (LNE)	Direct features: - on exposed soil, the occurrence of another spectral response (soil change) - texture change (cutting of existing textures on uncovered soil, dense vegetation) Indirect features: - other moisture of the former trench associated with the change in the soil structure (most often—more dryness than in the surroundings) - other vegetation/vegetation rate related to the violation of the original soil structure - geomorphological features (drainage network, dorsal network), the presence of paths, watercourses at the locations of the rampart/trench (provided that its layout was adapted to the existing terrain or the axis of the rampart/trench became a transportation route with time) - a contrasting spectral or textural border extending over a longer section between different areas	s (spectral) t (texture) m (moisture) v (vegetation) g (geomorphology), p (paths) w (watercourses) b (border)

.

Speaking about the Macedonian camp as a whole (encompassing an area of between 30 and 235 ha), it has also been assumed that the rectangular or square shape of the camp may have left remains of the whole structure or some parts of it, which would nowadays be visible as geometrized features (polygons). A detailed list of the photointerpretation features of the sought-after linear features with explanations is presented in

Table 3. Sought-after areal elements and related photointerpretation features anticipated in the satellite imagery.

Possible Current State	Photointerpretation Features, Associations, Indicators	Abbreviations in Figures
Agricultural activity, natural vegetation, local anomaly Geometrized remains of small features GIS encoding polygon vector P	Direct features: - the area inside the mapped contour of the camp rampart, with possible accompaniment or continuity through geomorphological elements of the terrain, such as paths or watercourses - shape, size—small-scale residues related to the inner layout of the camp Indirect features: - no contour, local anomaly characterized by a different spectral or textural response compared to the surrounding environment - no contour, local anomaly characterized by varied spectral response/texture, partly different from the surrounding environment	c (contour) sa (spectral anomaly) ta (texture anomaly) va (varied anomaly) paofr (prospective areas of research)

.

In analyzing the contents of Table 2 and Table 3, it is easy to notice that the features/anomalies were very diverse in nature, and the analysis required the use of extensive background material. Therefore, the time variation of recorded images (different seasons) and the simultaneous selection of multi-band and panchromatic recordings were adopted during the initial stage. This allowed for the use of appropriate, multi-variant image processing methods to enhance the image and facilitate interpretation:

• The development of different color composites (CC)
• Image data fusion/spatial enhancement: increasing the spatial resolution through the integration of multispectral and panchromatic images
• Use of principle component analysis (PCA) transformation
• Use of decorrelation stretch
• Use of vegetation indices

2.2.1. Development of Different Color Composites

The Pleiades satellites register 4 bands (B, G, R, NIR). In practice, two standard color composites are used for 4-band recording: CC in natural colors (B, G, R) and false color composite (FCC) in false colors, including a near-infrared band (G, R, NIR). These two composites were used in the interpretation process.

WorldView-2 records 8 bands. Combinations of triples of bands give a total of 56 possibilities, which, in turn, equals 336 color composites. For Beauchemin and Fung [103], the appropriate selection of bands for CC is based on the selection of spectral ranges that collectively contain the most information accompanied by low redundancy. To compare the CC and the degree of information contained therein, two statistical indicators were used, MOIK (named after its creator, J.G. Moik) (Equation (1)) (being the sum of the correlation modules of the three bands included in a given CC) [104] and the optimum index factor (OIF) indicator (Equation (2)) (showing the normalization of the sum of standard deviations and module coefficients, that is, absolute values/correlations of three pairs of spectral bands, from which a CC is formed) [105].

MOIK_CC(b,c,d) = |r_(b,c)| + |r_(b,d)| + |r_(c,d)|

(1)

OIF_CC(b,c,d) = (σ_(b) + σ_(c) +σ_(d))/(|r_(b,c)| + |r_(b,d)| + |r_(c,d)|)

(2)

Low values of the MOIK index indicate a slight redundancy of data for the given three bands. In turn, the OIF index is sensitive to the correlation of bands and additionally to the total amount of information contained in them—the higher the values of the indicator, the better. Neither of the two indicators solves the problem of arranging the three specific bands in the color composition (indicator values do not change), and the selection of the band order requires a visual analysis of the six available combinations each time. Both the band correlation parameters and the band variances for the analyzed WorldView-2 image are presented in

Table 4. Correlation of the WorldView-2 multispectral image bands.

Correlation	Band 1	Band 2	Band 3	Band 4	Band 5	Band 6	Band 7	Band 8
Band 1	1	0.9340	0.8825	0.8623	0.8356	0.7703	0.6717	0.6653
Band 2		1	0.9845 2	0.9715 2	0.9639 2	0.9128 2	0.8255	0.8236
Band 3			1	0.9916 1	0.9885 2	0.9620 2	0.8921	0.8902
Band 4				1	0.9937 1	0.9645 2	0.8825	0.8863
Band 5					1	0.9588 2	0.8814	0.8835
Band 6						1	0.9695 2	0.9719 2
Band 7							1	0.9939 1
Variance	497.3	2006.2	11136.2	22750.1	14316.8	25889.6	20833.5	16444.3

¹ Extremely high correlated bands; ² very highly correlated bands.

.

Rankings of composites according to the MOIK and OIF indexes, respectively, are given in

Table 5. List of information capacity of WorldView-2 (WV2) triplet of bands, ordered by: (a) optimum index factor (OIF) index and (b) MOIK (named after its creator, J.G. Moik) index.

(a) Position	Triplet of Bands	OIF	MOIK	(b) Position	Triplet of Bands	MOIK	OIF
1	467	161,927	2.81649	1	178	2.16059	136,477
2	468	155,867	2.82268	2	158	2.38432 1	113,316
3	478	153,252 1	2.76278	3	157	2.38862	119,853
4	567	151,225	2.80962	4	167	2.41145	135,825
5	457	150,427	2.75759	5	148	2.41391	124,845
6	456	147,887	2.91693	6	147	2.41655	131,372
7	678	147,676	2.93527	7	128	2.42287	80,617
8	367	145,477	2.82357	8	168	2.42944	128,191
9	568	145,260	2.81414	9	127	2.43113	86,966
10	347	144,850	2.76627	10	138	2.43791	105,033

¹ Selected triplet of bands (color composites).

. The results point to a large differentiation in the order of composites. This is due to the 11-bit recording of WorldView-2 (WV2) images, which translates into a very large variation of the standard deviation for individual bands (see Table 4). This results in the elimination of composites containing visible bands 1, 2, and 3 from the first places in the OIF index. In turn, color composites are preferred when there are very low correlations for two specific pairs of bands, while one of the pairs can be highly correlated in this case (e.g., CC 178 in the MOIK ranking; CC 478 and, 347 in the OIF ranking). These cases are undesirable, considering the aim of correctly selecting a CC. This means that these parameters, which work well for the Landsat data selection, malfunction in some cases for the WV2 dataset. Consequently, they should be treated only as an auxiliary tool in the selection of color composites, especially taking into account the critical duplication of information that is easy to observe in Table 4.

In analyzing the two indexes’ rankings (Table 5), correlations, and band variances (Table 4), and taking into account the visual qualities of the compositions, an attempt was made to find a compromise. The result of this compromise was that the color composites had relatively good places in the rankings and, at the same time, were based on bands with a relatively low correlation. In this case, both standard composites (used for the Landsat missions) and unusual bands present in WV2 (such as the coastal blue (band 1), yellow (band 4), red edge (band 6), and infrared NIR2 (band 8)) were taken into account in the selection procedure. Duplications of nearly/almost identical band pairs in the indicated composites were also avoided (for example, composite 178, occupying the first position in the MOIK index, would be similar to the selected composite 478 based on the MOIK index). Ultimately, a set of 4 color composites was used in further tests:

• CC235 (equivalent of TM123, natural colors, TC—true color)
• FCC357 (equivalent of TM234, CIR—color infrared: color composite: G + B + NIR)
• FCC478 (high positions in ranking the OIF, presence of unusual bands—yellow and red edge)
• FCC158 (high positions in ranking the MOIK)

2.2.2. Image Data Fusion/Spatial Enhancement

Two fusion methods (merging) were used to improve the spatial resolution of the WV2 image: high-pass filter (HPF) and intensity hue saturation (IHS) [106,107]. In the first step, each of the multispectral bands was subjected to preliminary processing, which involved adjusting to PAN resolution and consequently led to the removal of the blocking effect (see Figure 3).

The HPF method is characterized by the maintenance of the spectral characteristics of the input bands [108], and the degree of detail depends on the size of the high-pass filtration window. It involves the removal of thematic (spectral) information from the panchromatic image, leaving only details (edges) on it. This effect is achieved through HPF filtration, in which the sizes of the extracted structures and edges are controlled through the appropriate filter window. It is generally accepted that the filter window should be larger than twice the ratio of the combined images. In the next step, such an image is added to individual spectral bands (see Figure 4).

Because the results of HPF filtration are images with medium brightness equating to zero, this ensures that the original average brightness of the spectral bands is maintained in the HPF integration. After multivariate tests, the HPF filtration window of 9 × 9 was found to be the best (similar results were obtained for the 11 × 11 window). This is in accordance with the aforementioned theoretical principle of the selection of a filter window, since the ratio of the image resolution of MS to PAN equals 4:1.

For smaller windows, however, inconsistent spatial information was obtained with dominant large, primary pixels of the multispectral image, and the detail enhancement was too small. In the case of the tested larger filtration windows, spatial information dominated over the spectral—colored composites were deprived of color (reduced saturation), and the edges presented the main information in the picture.

The IHS method does not provide as much preservation of spectral information as the HPF method, but despite these reservations, it is characterized by high interpretative values. In this method, three selected spectral bands are enhanced spatially at once and are then used to build a specific color composite (see Figure 5). Therefore, the integration has been applied to specific triplets of bands, collated according to selected color composites: 158, 478, 235, and 357. Spectral images are subjected to transformation after pixel multiplication and after the removal of the blocking effect. Each time, after transformation from the RGB space to the IHS space, image I (intensity) was replaced with a PAN image, radiometrically corrected: its average brightness and variance were consistent with image I. This approach provided the maximum possible protection of the original spectral information of the bands, which could be distorted during the merging procedure.

The results of the HPF and IHS methods for color composite 158 (showing a particular case area) are summarized in Figure 6.

2.2.3. Use of Principle Component Analysis (PCA) Transformation

Principal component analysis (PCA) belongs to one of the procedures of image processing in which the creation of new images based on the input bands takes place through a weighted linear combination [76,78]. Standardized PCA was applied to all datasets, yielding four components for the Pleiades images and eight components for the WorldView-2 image. The basic tool for the statistical analysis of our data after PCA transformation was the loading matrix, i.e., the correlation matrix between the primary bands and the main components. It enabled an assessment of how the information flow took place from specific bands to the main components. Such information allowed for the appropriate selection of three components to obtain a diverse set of new color composites. On the basis of the analysis of the loading matrix, the color composites were developed based on the principal components:

• Pleiades: CC (PC1, PC2 PC3), CC (PC1, PC2, PC4)
• WorldView-2: CC (PC1, PC2, PC3), CC (PC1, PC4, PC5), CC (PC2, PC4, PC7), CC (PC6, PC7, PC8)

Based on the analysis of the loading matrix, attempts were also made to combine specific bands and main components in color composites. After visual tests, mixed composites, i.e., those containing components and selected bands, were disqualified from further processing (the mixed composites were rarely used and did not offer any important information; thus they were not included in

Table 6. Summary list of materials prepared for the visual interpretation.

Satellite	Pleiades	Pleiades	WorldView-2
Date of registration	11 April 2013	12.12+24.12.15	28 September 2014
PAN image with enhanced contrast (resolution 0.5 m)	No	No	Yes
Spectral bands with enhanced contrast black-and-white (BW) (resolution 2 m)	No	No	Yes, 8 bands
Spectral bands with enhanced contrast after merging (BW, resolution 0.5 m)	Yes, 4 bands	Yes, 4 bands	Yes, 8 bands (HPF and IHS)
CCs (RGB 24 bit; 0.5 m)	123, 234	123, 234	235, 357, 158, 478 (HPF and IHS)
Color composites with decorrelation stretch (RGB 24 bit; 0. 5 m)	123	123	235, 478 (HPF and IHS)
Principal components of the PCA (BW, resolution 0.5 m)	4 components	4 components	8 components (HPF and IHS)
Color composites containing principal components (RGB 24 bit; 0.5 m)	PC1, PC2, PC3 PC1, PC4, PC2	PC1, PC2, PC3 PC1, PC4, PC2	PC1, PC2, PC3; PC1, PC4, PC5 PC2, PC4, PC7; PC6, PC7, PC8 (HPF and IHS)
Vegetation indices Normalized Difference Vegetation Index (NDVI) and Perpendicular Vegetation Index (PVI) (8 bit)	Yes	Yes	Yes, (bands 5, 7 and 5, 8) (HPF and IHS)

).

2.2.4. Use of Decorrelation Stretch

Much better results were obtained when the method of decorrelation stretch was employed based on the PCA transformation [76,109]. This procedure is particularly useful for composites containing highly correlated spectral bands characterized by low contrast and poor coloring (e.g., a composite in true colors, equivalent to TM 123). This involves the standardization of the main components obtained after the PCA transformation, followed by the inverse transformation. The standardization of components results in exposing the unique spectral information located in further components and characterized by low variances. At the same time, information repeated in all bands in initial components with high variances is reduced. The effect is manifested through the enrichment and diversification of colors for features/areas that initially appear to be spectrally uniform. Decorrelation stretch was tested for various color composites. Finally, based on the activities described above, the enhancement for CC478 and CC235 was developed (see Figure 7), because for these compositions, the highest data redundancy was noted (presence of highly correlated bands)—see the MOIK ranking (see Table 5b).

2.2.5. Use of Vegetation Indices

In the present study, the Normalized Difference Vegetation Index (NDVI) and the Perpendicular Vegetation Index (PVI) were used. As was recently confirmed by Agapiou [94], good results were obtained from the use of both indices in satellite archaeology, especially the PVI. The NDVI belongs to indices from the group expressing the slope (proportions of brightness between specific bands), while the PVI belongs to indices expressing the distance from the determined soil curve. The NDVI is characterized by high dynamics (reaction to changes in the state of vegetation) and is resistant to so-called soil noise and atmospheric changes. The results are obtained on a linear scale in the range of −1 to 1 regardless of the radiometric range of the data. In turn, the PVI is dedicated to areas where, in addition to vegetation, a large proportion of soils (initial or final vegetation phase, semi-arid areas, etc.) are often recorded within the pixel. The index allows the reflection recorded from the uncovered soil and the reflection recorded from vegetation to be correctly distinguished, both located simultaneously on one pixel (mixel). At the initial stage, to calculate the index, it was necessary to determine, each time on each dataset, the soil curve by examining spectral reflection in locations of uncovered soil, characterized by varying degrees of moisture. Analysis of the NDVI and PVI images contributed to the detection of geometrized artifacts marked by vegetation patterns.

2.3. Interpretation of Imagery Data

After the initial processing, differentiated image products were obtained that were subjected to the process of photointerpretation. In total, for the three selected dates of recording, 57 black-and-white (BW) images were available (in grayscale), as well as 12 images in the vegetation indices scale and 30 RGB color images (see Table 6).

The interpretation process was carried out in ArcGIS (Environmental Systems Research Institute (ESRI), Redlands, California, United States). Individual black-and-white and color images were superimposed in the ArcMap module. Work was carried out mainly on color images, first using color composites (123 and 234 for Pleiades and 235, 357, 158, and 478 for WV2). Only later, due to identification difficulties, were other processed images used. The work was iterative: from general to detailed, with multiple scale changes and changes in background layers (use of different registration dates, exchange of observed images). A smaller-scale look at the area made it possible to capture the spatial context, linking the detected feature/phenomenon to others, especially the terrain. In some cases, this approach allows for the determination of the origin/function of the detected features or can ensure the existence of an anomaly that requires field verification. For large scales, the 0.5 m single-pixel size made it difficult to identify small structures. In such cases, it was useful to reach for true-color aerial photographs available in the Bing browser (for this area, Google Earth provides only VHR satellite data, i.e., with lower resolution than the Pleiades and WorldView-2 data). In many cases, the detected feature was visible only on specific materials (for example, it was visible thanks to plant cover, or vice versa, with the vegetation masking the feature). There was often uncertainty concerning the existence of the artifact in a given place—its presence was suggested by many images, but the identification was validated only by increasing the spatial resolution of data through the integration process (merging). Such situations confirmed the need for the acquisition of varied data and work on a complementary set of well-selected color composites. Unusual image processing, often difficult to interpret without resorting to classic color composites, also allowed for the identification of some local anomalies, usually associated with slightly larger areal inconsistent with land use. Such situations confirmed the desirability of the adopted methodology for preparing a diverse set of image data for photointerpretation (not only the simple creation of CC), including the need to have data from different seasons.

According to the anticipated, diversified photointerpretation features (see Table 2), features were usually mapped that showed several features indicating the possibility of the occurrence of the sought-after features. Embankments and trenches, in practice, were most often suggested by the presence of linear features (linear non-exist (LNE)) such as paths and landforms (areas of hills, visible on data from different seasons), linear changes in land desiccation (autumn and winter data), or increased vegetation (spring data). An additional indicator was texture disturbance, especially in the agricultural landscape, which facilitated the detection of linear features running through arable fields of various natures. Only occasionally were changes in the type of soil observed that would not be associated with changes in moisture content. Of all the mapped features of a linear nature, only in one case was a feature detected that could have features of a direct exposure of the rampart and trench (LNE) profile.

A tonal or color separation (less frequently textural) from the surroundings usually led to the mapping of surface features that could be traces of the existence of a camp (see Table 3). The observed separation could also be diversified within a given detection. The boundaries of this type of anomaly were fuzzy, often uncertain. In these cases, the precision of determining the border and shape of features was therefore of secondary importance.

Due to the low heights of features (such as buildings, embankments, trenches, and trees) and the relatively large image pixel (0.5 m), the role of shadows turned out to be irrelevant in the interpretation. Only on the hills and in the winter images could the role of shadows be slightly observed due to low illumination by the sun. The shadows highlighted the natural forms of the terrain, enabling mountain ridges and drainage networks to be more easily detected. The phenomenon of the area’s low illumination, however, turned out to be insufficient for detecting anthropogenic features, as is often the case when using aerial photographs in archeological prospection.

Many features were identified during the first stages of photointerpretation that required further verification. Therefore, in the process of the first verification, indirect features, mainly associations, gained importance and led to the elimination of features potentially related to the current land cover, existing transport and settlement network, plot boundaries, watercourses, landforms, etc. In many cases, deciphering connections between features (associations) allowed for the determination of their current character or function, and thus, they were excluded from further analyses.

The documentation of the photointerpretation analysis formed the basis for carrying out ground verification at a later stage. Thus, the approximate geometry of features detected in the images was sufficient in many cases. Situations of this type were frequent, especially for the category of areal features. As mentioned above, the features/anomalies in this category were characterized by inadequate outlines. Approximate polygons were also used for prospective areas of research so as to include small-scale features requiring further analysis. However, in the case of linear and areal features with distinct recognizable contours, vectorization was carried out precisely.

Documentation of photointerpretation works was carried out by vectorizing detected anomalies. The most popular shapefile vector data standard (ArcGIS, shp) was used. The adopted form of documentation allowed for the following:

• Self-correction was carried out on other sets of image data, which led to the confirmation of anomalies/features or the deciphering of their origin as not being related to the conducted research, and thus a reduction of identifications from the initial list.
• After positive verification, the documentation was sufficient to find a feature/prospective area in the field for final classification/verification.

The diversity of the nature of photointerpretation features of indications for the fieldwork was reflected in the vector documentation.

3. Results

Linear features can be identified relatively well despite their small width (even at a sub-pixel level). The continuity of this type of feature allows its course to be tracked, even despite the tonally or texturally changing surroundings. In the preliminary study area, 27 lines were determined after photointerpretation of a full set of image data (see Figure 8).

In hilly terrain, dorsal lines (lines connecting points along mountain ridges, and also lines separating adjacent drainage sub-basins) and lineaments (linear surface features, or the composition of the surface, which differ distinctly from the surroundings and presumably reflect subsurface phenomena [110]) constituted additional difficulty [111]. These could be misleading (as only natural features), but on the other hand, they could also be anthropogenic (as transportation routes or natural elements of a defensive or strategic character possibly related to boundaries of the sought-after camp; see Figure 9A). In the lowlands, an area covered by arable crops, the dominant feature for the winter and spring images allowing for the detection of linear features was a brighter phototone. Most often, such features were modern roads or remnants of pedestrian paths that, in unpaved terrain, were characterized by faster drying than the surroundings (arable fields; see Figure 9B). In a few cases, the detected linear feature was characterized by a darker phototone (which can be observed due to registration in the reflected infrared range). Most often, the reason for this is a change in water conditions (greater humidity of the area) and/or associated richer vegetation than in the surroundings (see Figure 9C).

A large part of the linear features visible in the images is justified in the current land management—its land cover, functions, geomorphological features, landforms, etc. As described in Section 3, attempts were made in the course of preliminary verification to eliminate them to limit the scale of subsequent ground verification. For example, in Figure 9A, arrows labeled A indicate geomorphological features of the area that may incorrectly suggest the existence of a trench/rampart. The local change in the brightness of discovered soils may have a diverse genesis, including a result of the existence of an old buried embankment/ditch (agricultural work may also bring to the surface different soil associated with over-drying that is invisible on the surface). Above all, the images reflect the current works of agricultural machinery, which may lead to misinterpretations (see arrows labeled B in Figure 9B). At the same time, in Figure 9C, arrows labeled C show an example of the presence of the foundations of a group of buildings with communication infrastructure, which may indicate the presence of old settlements in this area, possibly related to the space inside a camp. In fact, these are traces of recently abandoned modern buildings.

In accordance with the scheme adopted in Table 2 and Table 3 (interpretative features), the designated polygons have a dual character: they are to document the approximate course of detected areal features (areas) or anomalies as a potential camp area and indicate the places of existence of geometrized small-scale features that could be related to the internal infrastructure of the camp (see Figure 10).

Features/areas that had an areal character were designated if they could be distinguished from the surroundings by at least one photointerpretation feature. This group of features can be characterized by very large discrepancies in terms of size—from several dozen meters to even several hundred. In the area under study, after photointerpretation of the full set of image data, 30 polygons were finally determined. From these, 14 were located in the Mahad hills. Identified features can overlap each other, which results from their discovery in different images. Thus, the spatial range of a detected feature may have a completely different course than that of another adjacent feature detected in a different image; this highlighted other local diversity occurring in a given area.

The detected areal features were characterized by at least one photointerpretation feature that distinguished them from the surroundings. Their presence and intersection with the current state of land cover and land use resulted in complex projections, with their parts often disappearing over certain sections. In many cases, it was necessary to resort to a complementary use of various images from our database to determine potential artifacts. The following example shows the advantages of decorrelation stretch, based on PCA transformation (see Figure 11). Treating this type of product as support for a CC, in some cases one can see the distinctness of features that seem to be similar to others (see Figure 11C, with yellow and black arrows pointing to such an example). In turn, data generalization and the interpretation of tonal variability in a broader context offer an opportunity for detecting structures/features crossing traditional boundaries of types of land cover and land use or of visible plot boundaries (in Figure 11; white arrows indicate such a feature; the left lower limit of this feature was only tentatively estimated, while the remaining range can be distinguished by a relatively clear contour. The anomaly was intentionally not marked as a polygon to visualize the interpretation problems encountered in determining the course of fuzzy boundaries).

4. Discussion

The focus of this study was on the discovery of the possible remains of the Macedonian military marching camp located on the hills west of the modern town of Mahad (see Figure 2). The Macedonian marching camp was a space enclosed by a protective circuit trench and a low wall of dug soil thrown up behind the trench. Given the fact that the camp was used by 47,000 Macedonian troops, it must have enclosed a space of considerable size, approximately between 30 and 192 ha (different estimates).

Indeed, the question of the size of the sought-after Macedonian camp turned out to be most decisive in our research. Of the 14 polygons identified on the Mahad hills, the size of most of the structures did not even come close to the expected sizes. The two biggest polygon structures measured only 15.2 and 25.4 ha, respectively (see Figure 12). Theoretically, one could consider whether a structure of such size could actually have been part of the inner layout of the camp (for instance, the royal headquarters center cordoned off the main camp—[52]). However, due to the lack of detailed historical or archaeological data on Macedonian camps, one would remain in the sphere of pure speculation when suggesting such structures.

Apart from the aforementioned surface features, the designated linear features were an important element in the search for structures or their fragments. Nineteen such objects have been identified within the Mahad hills. Their number only allows for preliminary conclusions as to the usefulness of the collected image materials in the search for this type of field anomaly. Despite the above-mentioned reservations, the obtained results (19 detected linear features) allowed us to critically evaluate our materials and methods, as well as image acquisition dates. Consequently, we were able to create the appropriate rankings (

Table 7. Ranking of usefulness of multispectral data according to the acquisition times (seasons).

Ranking	Season	Unique Weight = 4	Primary Weight = 2	Secondary Weight = 1	Score ∑ = U + P + S
1	IV (autumn)	4	6	4	32
2	II (spring)	2	2	9	21
3	I (winter)	0	5	4	14

,

Table 8. Ranking of usefulness of photointerpretation features in the analysis of multispectral data.

Ranking	Photointerpretation Feature According to Table 2	Unique Weight = 4	Primary Weight = 2	Secondary Weight = 1	Score ∑ = U+P+S
1	Moisture	1	5	7	21
2	Spectral	1	4	1	13
3	Texture	0	1	8	10
4	Geomorphology	0	3	2	8
5	Paths	0	2	4	8
6	Vegetation	0	1	2	4
7	Border	0	1	1	3

and

Table 9. Ranking of the usefulness of the multispectral data processing methods used by the authors.

Ranking	Materials	Primary Weight = 2	Secondary Weight = 1	Score ∑ = U + P + S
1	CC color infrared (CIR)	6	5	17
2	CC WV2 158	2	12	16
3	CC Pleiades PC1, PC2, PC3	3	10	16
4	CC WV2 478	2	10	14
5	CC WV2 PC1, PC2, PC3	2	9	13
6	DS WV-2 478	2	6	10
7	DS TC	1	7	9
8	CC TC	1	6	8
9	CC Pleiades PC1, PC4, PC2	1	5	7
10	Vegetation indices	0	5	5
11	CC WV2 PC1, PC4, PC5	0	3	3
12	CC WV2 PC2, PC4, PC5	0	3	3
13	CC WV2 PC6, PC7, PC8	0	2	2

). These rankings may offer practical tips for future work in studies of this type, especially in the initial stages of selecting imaging materials.

Table 7 presents a ranking of the data in terms of their acquisition in a specific season of the year. As the detected linear features could be identified on one or several datasets, the evaluation adopted the weighting method in calculating the score. In unique cases when the detection of the features was possible only at a specific time of the year, the highest weight was assigned (of unique importance translating into weighted value = 4). When it was possible to identify the features in images from several seasons, the scores were differentiated depending on when it was easiest to identify the linear feature (of primary importance translating into weighted value = 2) and whether there were secondary indications confirming the existence of the feature (of secondary importance translating into weighted value = 1). The sum of the points determined the order of usefulness of the data.

A similar ranking was constructed (see Table 8) according to the importance of photointerpretation features described in Table 2. The results of the photointerpretation in graphic form were presented in Figure 12.

The third ranking concerned the usefulness of the multispectral data processing method used to detect 19 linear objects (Table 9). Only selected groups of materials were included in the ranking list; those excluded were difficult to interpret or generally useless (the following materials were excluded: individual spectral bands, color composites in the original 2m resolution, PAN image, components of PCA transformation). No distinction was made between the data after the IHS and HPF merging because, in practice, these materials yielded very similar results. The column about unique identifications was not included in the ranking as such a case never occurred (linear objects detected in the photointerpretation process always appeared on more than one set of image data).

The suggested features were subjected to field verification between 3 and 15 September 2019. The field groundtruthing did not observe any diagnostic signs of anthropic features (anthropic geomorphological features/abnormalities in the original terrain that cannot be attributed to orogenetic processes) in the terrain or concentrations of archaeological finds and, consequently, ended in negative results. In fact, the areal features were concluded to be natural rock formations. Likewise, of the 19 linear features identified on the Mahad hills, most turned out to be traces of modern transportation lines or natural phenomena. Furthermore, they did not form any reasonable perimeter characteristic of enclosures. Only in one case did the present survey come across traces of undoubtedly ancient remains—the line identified east of the Mahad hills (see Figure 12 and Figure 15), running roughly north-south, accounts for a short stretch of the over-90-km-long Assyrian canal system built by the Assyrian king Sennacherib (704-681 BCE) to bring water from Khinis at the footsteps of the Kurdish mountains to Nineveh [44,112,113,114,115,116,117].

The determining factor leading to the detection of the Assyrian channel in our images was the diversity of data acquisition in terms of seasons. As can be seen in Figure 13A–C, the channel is visible mainly on the autumn images. In detecting the course of the canal in our images, both vegetation and local changes in water conditions—local land desiccation—proved to be important indicators. These factors became important mainly in the autumn season, when the persistence of stable weather led plant vegetation into its final stage and agricultural work was completed (harvest in most fields). As a result, a linear feature became visible in the area of uncovered soils, thanks to its thick plant cover (black arrows in Figure 13C). In other nearby locations, its linear course was also reflected through a clearer spectral response resulting from a lower level of soil moisture (white arrows in Figure 13C). By analyzing all these detected places, it became possible to fragmentarily map the course of one linear feature.

In the case of the data from the other two vegetation periods, either the linear feature was masked by intensive vegetation (spring registration, Figure 13B), or vegetation disappeared with a simultaneous leveling of local differences in humidity in the area of uncovered soils (winter registration, Figure 13A). On both of these datasets, the channel course can be mapped (e.g., in the areas indicated by the arrows in Figure 13A,B), but the information is too meager to reliably determine that it is one continuous linear feature of considerable length. However, the winter and spring material helped complement the analysis of the autumn images, authenticating a local interpretation of images and confirming the detection of the feature. In addition, it should be noted that the enhancement of the images (decorrelation stretch and vegetation indices; see Figure 13D,F) did not result in valuable interpretive material for the examined feature (the channel was partly visible only on the color composites developed based on the principal components; see Figure 13E).

It should, however, be noted that in this part of its course, the canal is covered underground and is recognizable on the surface, if at all, only through indirect photointerpretation features (see Figure 14A), unlike in other parts of its course, for instance at Jerwan (see Figure 14B).

In this respect, the question has been posed as to whether the Assyrian canal, which encircles the entire eastern and southern border of the hills in question, could have been used as a trench for the purpose of setting up the Macedonian camp (see Figure 15).

In fact, it has been suggested that the Mesopotamian water canals could also have been occasionally used for military purposes as barriers [117]. However, given the fact that the Assyrian canal is placed below the Mahad hills, it would have meant placing a military barrier below the uphill positions taken by the Macedonians. While not impossible, this would have been an example of very irregular tactics—one would normally expect embankments to be placed uphill in such terrain to gain more height and distance from potentially approaching enemies.

5. Conclusions

Summing up the results of our research focused on the photointerpretation of VHRS data so far, it can clearly be seen that the independent potential of the data is insufficient to solve the problem at hand. Despite the use of various techniques of multispectral image enhancement, which can now be carried out on a digital form of data (i.e., broader than, for example, those used in CORONA data analysis), the vast majority of suggestions made within the course of remote sensing research were not confirmed during fieldwork. As for the photointerpretation, it would be advisable to strengthen the analysis by using detailed DEM models and SAR imaging.

Despite the failure to achieve the main goal, it is possible to provide some initial guidance on the use of multispectral data as an essential element for this type of research. The analyses confirmed that the best interpretative material is the autumn data, followed by the spring data (see Table 7). The most convenient interpretation work was carried out on color compositions, while the products of image processing (PCA, decorrelation stretch, NDVI, PVI) were additional material (see Table 9). An important issue is to enhance the spatial resolution from 2m to 0.5m by fusing the images (materials in the original resolution of 2m turned out to be useless). The choice of integration method was, in the end, irrelevant (comparable results were obtained for both IHS and HPF methods). Among the tested color composites, the highest suitability for CIR compositions was observed, which confirmed the results in similar satellite archaeology studies (focusing on R and NIR registrations) published by other scholars [74]. An almost identical result was obtained for CC WV-2 158, and a slightly lower result was obtained for CC WV-2 478. This indirectly confirms that the selection of composites for WorldView-2 according to the formal indicators OIF and MOIK was chosen correctly by the authors. In practical terms, however, this means that the additional bands of the WorldView-2 satellite do not offer any additional benefits when compared to the four-band recording of Pleiades (RGB + NIR). In contrast, the use of PCA transformation, as well as the decorrelation stretch used on visible bands (the third and seventh ranking positions, respectively, according to Table 9), may offer considerable benefits and may consequently be used to supplement the use of CC CIR. Due to the time and cost of data processing, it is recommended to use the multispectral Pleiades image with a resolution of 0.5m for further research. The most important photointerpretation feature was the local differentiation of humidity, followed by the spectral and textural characteristics (see Table 8). Therefore, it is also worth considering the use of times-series data in the future, which may highlight humidity anomalies occurring periodically or even episodically. Our observations are in line with the conclusions regarding the construction of a hypothetical optimal satellite system for archaeology [118].

In terms of archaeological research, this study reports the first results of our search for the material remains of the Gaugamela battlefield, specifically the remains of the Macedonian military camp set up on the hills directly overlooking the battlefield on the eve of battle. It is concluded that despite the use of very high resolution multispectral images, no traces have been found that could be interpreted as remnants of an earthen enclosure capable of accommodating around 47,000 soldiers. Several possible explanations and future recommendations can be given. First, it is theoretically possible to say that the negative results speak against the identification of the Navkur Plain as the location of the Battle of Gaugamela. At the same time, it should be noted that no such traces have been recorded for alternative locations either, especially in the case of the Karamleis–Qaraqosh hypothesis. Second, it is also possible that the wrong area in the Navkur Plain has been chosen, and, for instance, the search area should be enlarged to include the area north of Tell Betnar and southeast of Tell Mousakan (see Figure 16). While both areas are at a lower altitude than the stretch of hills between Tell Betnar and Tell Mousakan, they are still at a considerably higher altitude than the plain west of Tell Gomel (Gir-i Gomel).

Third, it is also likely that there are a number of reasons why the remains of the Macedonian marching camp have not been preserved. In this context, it should be noted that the Macedonian camp was a marching camp and not, for instance, a winter camp. A winter camp would have had a longer occupation, which may, in turn, have left more durable traces. In this context, the nature of the military enclosure at Gaugamela must be contrasted with the nature of military enclosures in Dura Europos and Hatra. Namely, in both cases, the findings were remains of more massive structures that were used for a much longer period of time. Hatra’s circumvallation walls (including both stone and earthwork walls and ditches) were most likely created by the Sasanian army, which besieged Hatra for several months in 240/241 CE [119]. In turn, three enclosures at Dura Europos (both brick and earthwork, no ditches) may have been in use on a number of occasions between 115 and 256 CE [67]. Finally, the Navkur Plain is an alluvial plain where river deposits may effectively remove or mask archaeological evidence. The Navkur Plain has also been extensively used for agriculture by local populations, unlike the desert areas of Hatra and Dura Europos, and traces of ancient entrenchments may have been eradicated by centuries-long agricultural activities [44,120].