Detecting Classic Maya Settlements with Lidar-Derived Relief Visualizations

Abstract

In the past decade, Light Detection and Ranging (lidar) has fundamentally changed our ability to remotely detect archaeological features and deepen our understanding of past human-environment interactions, settlement systems, agricultural practices, and monumental constructions. Across archaeological contexts, lidar relief visualization techniques test how local environments impact archaeological prospection. This study used a 132 km² lidar dataset to assess three relief visualization techniques—sky-view factor (SVF), topographic position index (TPI), and simple local relief model (SLRM)—and object-based image analysis (OBIA) on a slope model for the non-automated visual detection of small hinterland Classic (250–800 CE) Maya settlements near the polities of Uxbenká and Ix Kuku’il in Southern Belize. Pedestrian survey in the study area identified 315 plazuelas across a 35 km² area; the remaining 90 km² in the lidar dataset is yet to be surveyed. The previously surveyed plazuelas were compared to the plazuelas visually identified on the TPI and SLRM. In total, an additional 563 new possible plazuelas were visually identified across the lidar dataset, using TPI and SLRM. Larger plazuelas, and especially plazuelas located in disturbed environments, are often more likely to be detected in a visual assessment of the TPI and SLRM. These findings emphasize the extent and density of Classic Maya settlements and highlight the continued need for pedestrian survey to ground-truth remotely identified archaeological features and the impact of modern anthropogenic behaviors for archaeological prospection. Remote sensing and lidar have deepened our understanding of past human settlement systems and low-density urbanism, processes that we experience today as humans residing in modern cities.

1. Introduction

For nearly a century, remote sensing has been used for archaeological prospection to locate otherwise hidden features and better understand past human-landscape interactions. While advances in higher quality and more easily accessible remote-sensing data have drastically increased our ability to conduct archaeological prospection [1,2,3,4,5], they also highlight how landscapes, vegetation, and visualization tools directly impact our ability to detect small archaeological features [6,7,8,9,10]. Particularly in the past two decades, novel technologies such as Light Detection and Ranging (lidar) revolutionized the use of remote sensing for archaeological prospection. In the tropics, these technologies can improve archaeological prospection through the remote and digital identification of archaeological features [11,12,13,14,15,16], thus increasing the efficiency of pedestrian archaeological survey [17,18]. Here, I assess the applicability of four methods (three relief visualization techniques and an object-based image analysis (OBIA)) based on aerial lidar-derived data, to remotely detect small archaeological features associated with hinterland settlements in the Maya Lowlands.

Archaeological features must be “reasonable representations of on-the-ground features” for archaeologists to use remotely sensed features in analyses [19] (p. 7). While lithic scatters and other ephemeral archaeological sites may be nearly impossible to identify by using remote sensing, the remains of small and large archaeological and historic masonry or earthen features can be recognized by using remote sensing [20,21,22,23]. Specifically, lidar inherently changed the nature of archaeological survey in tropical regions [11,12] and is used to identify archaeological features on a variety of geographic and environmental landscapes. Remote-sensing methods have been used to test lidar-derived models for archaeological prospection in a variety of spatial and temporal regions, including the Archaic period (pre-1000 BCE), Woodland Period (1000 BCE–1000 CE), and Mississippian Era (900–1200 CE) of the Eastern and Southeastern United States [21,24,25,26,27,28]; Bronze Age (2000–800 BCE), Iron Age (800–1 BCE), and Medieval (800–1300 CE) Europe [29,30]; Preclassic (1000 BCE–250 CE) and Classic period (250–800 CE) Maya in Central America [1,2]; Khmer Empire (800–1400 CE) in Cambodia [12,31]; and the Four Corners region of the American Southwest (600–1300 CE) [32].

The application of aerial lidar in archaeological contexts has gained immense popularity in the last decade and has fundamentally shifted the effectiveness of archaeological survey [11] (see Reference [33] for an in-depth review of the use of lidar within archaeological contexts). Archaeologists use a variety of spatial analyses and relief visualization techniques to aid in archaeological prospection. These techniques include hillshades of varying illumination angles [6,9,34,35], slope models [6,9,19,35,36], local relief models and simple local relief models (LRM/SLRM) [1,35,36,37,38], sky-view factor (SVF) [1,6,17,29,30,35,36,39], topographic position index (TPI) [40,41,42], geomorphon [33], topographic openness [29], prismatic openness [1], Principal Components Analysis (PCA) [6,17], Red Relief Image Maps (RRIM) [1,8,25,43], and OBIA [8,21,25,44].

Lidar analyses are increasingly popular [45], especially in areas with dense vegetation such as the Maya Lowlands, where traditional satellite imagery produced challenges for archaeological prospection [2,46]. Lidar produces a high-resolution Digital Terrain Model (DTM) of the bare earth, allowing users to detect minor features on the landscape, without the hindrance of vegetation. The quality of lidar data and ability to detect small archaeological features vary depending on the density and height of vegetation, as a consequence of modern human activities [7,8,9,10,41], and will affect the resulting topographic visualizations and computer-based learning analyses [21]. In the tropical environments of the Maya Lowlands, smaller archaeological remains, such as house mounds, have proven more difficult to detect through satellite imagery [47] and lidar [9,17,41,42]. Here, I compare relief visualization techniques for the identification of archaeological features and how the size of the feature, vegetation height, and modern land use impact remote archaeological prospection.

The present study focused on the visual identification of small, Classic (250–800 CE) Maya settlements, dispersed across an upland landscape in Southern Belize (Figure 1). Between 2005 and 2018, pedestrian settlement survey was conducted at two Classic Maya centers, Uxbenká and Ix Kuku’il, covering approximately 35 km² [48]. This survey identified more than 230 household groups [49,50,51] that are further divided into 339 discrete plazuelas. Plazuelas contain archaeological features such as house mound structures and public architecture situated around a central plaza. In 2011, the Uxbenká Archaeological Project (UAP) obtained a 132 km² swath of lidar coverage [9]. Pedestrian archaeological survey is expensive and time-consuming [42] (pp. 351–353), and only a fraction (26.5%) of the UAP lidar area has been surveyed. An extensive pedestrian survey is yet to be conducted between the major Maya centers of Southern Belize, and this study informs the extent and density of previously undocumented hinterland settlements, specifically between Uxbenká and Lubaantun (Figure 1). Here, I expand on previous studies [9,18,52], to compare three topographic relief visualization techniques—SVF, TPI, and SLRM—and OBIA on a slope model, derived from a high-resolution (1 m DEM) lidar dataset for the detection of Classic Maya plazuelas and archaeological features associated with hinterland household groups. SVF, TPI, and SLRM relief visualization techniques were selected for this analysis because they proved successful in the identification of Maya house mounds and plazuelas among other Classic Maya centers, using lidar datasets [1,8,41]. OBIA was selected because, recently, Davis and colleagues [25] successfully used OBIA and segmentation on lidar-derived relief models to identify mounded archaeological features. While OBIA is traditionally applied to satellite imagery, using OBIA on a slope model is a novel application that allowed me to circumvent the dense vegetation and apply the concept of segmentation to identify archaeological features (i.e., plazuelas) on the lidar data. The two techniques where archaeological features were easiest to visually detect on the UAP lidar (i.e., not using automated feature extraction) were TPI and SLRM. These two techniques were used to identify archaeological features (Classic Maya plazuelas) and then compared to archaeological features previously documented during pedestrian survey. The impact of vegetation height and modern land use on the ability to remotely detect archaeological features is also evaluated.

These findings suggest hundreds of new possible plazuelas exist within the 132 km² lidar zone. The results show variations in the SLRM and TPI for the non-automated visual identification of small archaeological features. When compared to the pedestrian survey data, the remotely identified plazuelas had a false-positive rate of approximately 14.7%. The results of this study add to the discussion of the utility and drawbacks of using remote sensing for archaeological prospection. Unsurprisingly, the findings of this study suggest that larger plazuelas and archaeological features are generally more discernable in the lidar-derived relief visualizations, although this relationship is inconsistent among smaller features, due to the size of structures and plazuelas and vegetation cover. This study emphasizes both the utility of relief visualization for the detection of archaeological features and the importance of ground-truthing and pedestrian survey [17]. These findings further our understanding of the extent and density of Classic Maya settlements in previously undocumented survey zones, highlighting the variations in the relief visualization techniques, continued need for pedestrian survey in conjunction with remote sensing, and the utility and shortcomings of remote archaeological prospection.

2. Materials and Methods

2.1. Study Area

Uxbenká and Ix Kuku’il are in the southern foothills of the Maya Mountains of Belize. The Southern Belize foothills were home to several large polities including Uxbenká, Ix Kuku’il, Lubaantun, Xnaheb, Nim Li Punit, and Tzimin Che, some of which erected carved and dated stone monuments [53,54,55] (Figure 1). The topography of Southern Belize varies from mountain and foothill uplands, with elevations between 250 and 500 m above sea level, to flat coastal plains. The topography largely influenced the location of Classic Maya settlements, as most sites are situated on areas of hilly relief in the foothills of the Toledo Uplands [56], which comprise only 6% of the Southern Belize landscape [57] (p.169).

2.2. Dataset

The National Center for Airborne Laser Mapping (NCALM) at the University of Houston collected aerial lidar data in May 2011 (see Reference [9] for details on the UAP lidar acquisition). The UAP dataset averages 20.1 returns per m² (for all returns) for the core of the dataset, and 13.8 returns per m² on the periphery of the dataset. The UAP data have an average of 2.72 ground returns per m² [9]. Ground returns for other lidar datasets in the Maya region are similar, ranging from 1.1 to 5.3 ground returns per m² (gr/m²) in Central Petén [1], to 1.35 gr/m² at Caracol [58], to 1.5 gr/m² at Yaxnohcah [59], to 2.8 gr/m² in the Belize Valley [60] (p. 8674), to 2.84 to 42.84 gr/m² at Ceibal, based on varying vegetation height [43]. The differences in ground returns are due to variations in sensors [46,61] and vegetation cover and density. Variations in vegetation include semi-urban and developed areas [41], shifting agricultural cycles [6,8,9], vegetation types in different microenvironments [10], and heavily forested regions such as preserves, reserves, and parks that have not undergone major anthropogenic modifications for decades [2,8,17]. Old forest growth on protected lands results in cleaner point clouds that facilitate the ability to detect archaeological remains on the resulting relief visualizations [62]. The number of ground returns directly impacts the ability to detect archaeological features. In general, higher ground returns (5–10 gr/m²) result in the detection of more archaeological features, while lower ground returns (1 gr/m²) yield lower detection rates [63]. Combined, the relatively low ground returns (2.72 gr/m²) and variations in vegetation height and density in the UAP lidar dataset will likely decrease the ability to detect some archaeological features. However, these variations also enable the assessment of how vegetation height and land use impact archaeological prospection on different relief visualizations.

2.3. Previously Detected Archaeological Features

The UAP conducted pedestrian archaeological settlement survey at Uxbenká and Ix Kuku’il from 2005 to 2018 [50,51]. The efficiency of the survey drastically increased with the acquisition of lidar data in 2011 and the ability to target minor hilltops for survey [9,18]. The survey identified a total of 339 plazuelas, which are spatially discrete groupings of archaeological features, including house mounds or monumental architecture, situated around a central plaza and often on an isolated hilltop or knoll of a hill (Figure 2). Of these, 315 were located within the lidar dataset. The survey identified an additional 68 hilltops in the lidar area that lacked archaeological features, which provided a basis for assessing the visual appearance of hilltops without archaeological features. The UAP survey data were used to test the effectiveness of remotely sensed archaeological prospection for both confirmed positives where known plazuelas are located and false positives where pedestrian survey identified no archaeological features. Most of the previously documented archaeological features are composed of ancient house mounds rather than monumental architecture.

Previous studies using the UAP lidar data [9,18,52] focused on hillshades, slope models, and LAS point clouds to identify Classic Maya house mound structures, concluding that individual structures are difficult to detect using a hillshade. In the 4 km² area surrounding the Uxbenká site core, only three of the 135 documented structures were visible on the hillshade [9] (p. 8). However, modified hilltops are more easily detected with a slope model: 56.5% (n=13) of 23 documented settlement groups near the Uxbenká site core were visible on the slope model [9]. The previous analyses focused on a 4 km² subsample constituting approximately 3% of the lidar data; this study is unique in that it expands on the previous research to areas of hinterland settlements across the 132 km² area, including lands that have yet to be surveyed, and appear diverse in regard to plazuela and architectural size and complexity, similar to the settlement systems of Uxbenká and Ix Kuku’il.

2.4. Relief Visualization Techniques

Different relief visualization techniques work better in certain contexts, based on local topography and vegetation, which impact the lidar point cloud and resulting raster files [6,8,41,64]. I used three relief visualization techniques (SVF, TPI, and SLRM) plus OBIA (on a slope model) on the UAP lidar dataset, with the goal of identifying new archaeological features. The average Classic Maya house mound in Southern Belize is small, varying from 20 to 275 m², and is less than 1.2 m in height [65], although the majority of house mounds are only a few courses of stone tall and generally less than 0.5 m in height. Non-domestic structural platforms are larger, with greater basal areas and heights. The tallest documented building at Uxbenká is 12 m high [9] (p. 3). The average household plazuela group is between 20 and 40 m in diameter, while the civic ceremonial plazuelas are upwards of 100 m in diameter.

To test for the ideal parameters in the identification of archaeological features, I ran several iterations of each tool with different inputs (see Figure 3 for final inputs). After the raster images were created from the iterations of analyses, I selected the best raster image from each method for identifying archaeological sites. Like others (see Reference [64]), to improve the visual detection of archaeological features, I modified the color ramp display for each raster output. Then, I visually assessed each raster image for its utility to identify Maya plazuelas and house mounds. I visually examined un-surveyed areas and previously surveyed plazuelas, to determine which outputs were best for identifying new plazuelas.

Of the four methods, SLRM and TPI produced results conducive to non-automated visual identification of modified hilltops and structures of the Classic Maya. Using a red (high) to blue (low) stretched color ramp classification for both the SLRM and TPI, I identified plazuelas based on the presence of flattened and modified hilltops which generally appear as a yellow plaza surrounded by a ring of red (Figure 4e). Plazuelas were also identified by the appearance of several red structures situated around a central plaza, even if the hilltop modification was less visible. Possible looting activity, indicated by small “donuts” of a dark red outline with a lighter red or yellow hole in center (Figure 4f), also suggested the presence of archaeological features.

After I visually identified modified hilltops on the SLRM and TPI raster outputs, I compared the newly identified plazuelas in each method (e.g., SLRM) to the following: (1) the other method (e.g., TPI); (2) previously surveyed and documented plazuelas (residential and public groups); (3) areas that contained no archaeological features; and (4) the impacts of modern anthropogenic activities (i.e., varying height of vegetation because of swidden or milpa agriculture) on the ability to remotely detect archaeological features (Figure 3). No additional analyses were conducted on the OBIA and SVF raster outputs.

2.4.1. Object-Based Image Analysis (OBIA)

The application of OBIA to satellite imagery, and more recently lidar data, allows for the detection of object-based features based on user inputs [66]. OBIA uses segmentation to distinguish individual objects based on grouped pixels or features. The strength of OBIA is its ability to detect and delineate geographic features, including buildings, streets, and vegetation. In recent years, the detection of objects and features, using OBIA on lidar-derived data, has proved successful in archaeological contexts (see References [3,67]). While OBIA was not developed as a relief visualization technique, it can be used on relief models. Specifically, OBIA was conducted on lidar data, to identify archaeological features in the Southeastern United States [21,25].

Here, I applied the concepts of OBIA and segmentation to detect and delineate flattened and modified hilltops and structures on a lidar-derived slope (degrees) model. On the slope model, the flattened, modified hilltops are visibly different than the surrounding hilly landscape [9], allowing for the delineation of these anthropogenically modified features. OBIA can be conducted in several software programs, including eCognition, GRASS GIS, and Esri ArcMap [66]. Using the Mean Segmentation tool in ArcMap, I set the spectral and spatial details to 15 and the minimum segmentation size to 20 pixels (m). The output raster delineated flat spaces, including modified hilltops and some structures (Figure 4d), but also delineated other flat spaces, such as valley bottoms. Other applications of OBIA to the UAP lidar dataset may yield more productive results.

2.4.2. Sky-View Factor (SVF)

SVF is a visualization technique that overcomes sunlight illumination issues present in typical hillshades [33,39,68,69] and provides the relative elevation of neighboring points, allowing the user to identify small relief features [30]. This technique is useful in archaeological contexts because it identifies steeper anomalies and corners associated with structures [35] (p. 238). SVF was calculated, using the open source Relief Visualization Toolbox (RVT [70]), version 1.3 [30,68,71]. Inputs for the SVF model included a GEOTIFF of the Digital Elevation Model (DEM), the number of directions, the search radius, and options for noise removal. Other archaeological studies used 16 directions and a 10 m search radius to identify small-to-medium anthropogenic landscape modifications in Slovakia [29], and the same inputs were suggested for larger hillforts in Slovenia [30]. Kokalj and colleagues [71] suggest using a search radius between 5 and 10 pixels (cells/meters) for smaller archaeological features. Because the Classic Maya plazuelas and house mounds in Southern Belize are much smaller than Roman Hillforts, I ran several iterations of the tool with various inputs, but ultimately used the raster output with 16 directions (the optimal number of directions according to [71]) and a search radius of 5 m.

2.4.3. Topographic Position Index (TPI)

TPI reflects the difference in values of a central pixel and the mean of the pixels around it [72] and relief variations in both large- and small-scale landscapes, based on the neighborhood input parameters [41]. TPI is calculated manually, using the Raster Calculator in ArcMap [73], or using an open-source toolkit (Topography Toolkit [72,74]). GIS4Geomorphology [73] describes the raster calculator method where the user creates three new DEMs (minimum DEM, maximum DEM, and smoothed DEM) based on a desired neighborhood size, using the Zonal Statistics tool within Spatial Analyst [73]. After the three new DEMs are created, a simple raster calculation is conducted, using the following formula [73]:

(“10 × 10” − “minDEM”)/(“maxDEM” − “minDEM”)

where 10 × 10 is the size of the cells used in focal statistics.

The final TPI raster is then manipulated to create the ideal imagery based on color ramp shading and classifications for the given dataset. In my analysis, this method, the raster calculator method, did not work well on the hilly topography of Southern Belize, but perhaps produces better results on datasets with less topographic variations, such as plains and floodplains.

The second TPI method uses the Topography Toolkit [72,74], which is freely available to download. I used several inputs for the search radii of the height and weight (e.g., 33 and 33; 10 and 10; and 2 and 10) to test which variables resulted in the best TPI image for archaeological features including plazuelas and structures. The Topography Toolkit TPI method produced excellent results for the UAP lidar dataset.

2.4.4. Local Relief Model (LRM) and Simple Local Relief Model (SLRM)

LRM/SLRM identifies local variations in topography, essentially flattening extreme variations on the landscape such as peaks and valleys and highlighting minor landscape variations [33] (p. 93) and is ideal for small variations in elevations, such as building foundations, structures, and depressions [1]. There are several methods to produce an LRM/SLRM dataset, including a free LRM toolbox download [75], a model builder work-flow [38], and the RVT version 1.3 [30,68,71]. Like other archaeologists in the Maya region [1], I used the RVT to calculate the SLRM. Kokalj and colleagues [71] suggest using an input radius that is slightly larger than half of the size of the feature that one is trying to identify. I used “10” for the radius for trend assessment in 1-x-1-m pixels because the goal was to identify Classic Maya house mounds and small hinterland plazuelas (Figure 3).

3. Results

The visibility of Classic Maya plazuelas and archaeological features varied in the four methods. While other studies successfully used OBIA [21,25] and SVF [1,35] to identify archaeological features, the OBIA and SVF produced results where plazuelas and structures were difficult to identify on the UAP data (Figure 4c,d). Therefore, no further analyses were conducted on the OBIA or SVF datasets.

3.1. Topographic Position Index (TPI)

The TPI raster was calculated using two different methods. I used the previously described Raster Calculator method [73] with 10 × 10 m cells to calculate the TPI. The outcomes of this method did not produce a raster image where archaeological features were easily visible (Figure 5c). The second method used the TPI function within the Topography Toolkit [72,74]. I inputted several heights and widths (Figure 5d–f), to produce the best visual image. The best TPI output (parameters: neighborhood: rectangle; height: 33; width: 33) clearly shows flattened, modified hilltops and structures (Figure 5f), and therefore was used to visually identify Classic Maya plazuelas.

In total, I identified 503 possible plazuelas using TPI across the UAP lidar zone. Of the 503 possible plazuelas, 381 were also identified by using SLRM (see below), and 122 were only identified on TPI (Figure 6;

Table 1. Comparison of plazuelas identified with SLRM and TPI.

	All TPI	All SLRM	Total all SLRM and TPI
Total plazuelas identified	503	580	702
Number of confirmed positives (Previously surveyed plazuelas identified)	111	117	129
Number of the 315 surveyed plazuelas not identified	204	198	186
Percent of 315 surveyed plazuelas identified	35.20%	37.14%	40.95%
Total newly identified plazuelas	385	457	563
Number of false positives (Previously documented hilltops with no archaeological features)	7	6	10
Percent of 68 previously documented hilltops with no archaeological features	10.30%	8.80%	14.70%

and

Table 2. Detailed comparisons of results from only SLRM, only TPI, and both SLRM and TPI, as well as all plazuelas identified with SLRM and TPI (see Table 1).

	SLRM Only	SLRM and TPI	TPI Only	All SLRM	All TPI	Total of all SLRM and TPI
Plazuelas identified	199	381	122	580	503	702
Confirmed positives	18	99	12	117	111	129
Newly identified plazuelas	178	279	106	457	385	563
False positives	3	3	4	6	7	10
2011 milpas	14	39	1	53	40	54
2010 milpas	8	12	0	20	12	20
2009 milpas	6	9	2	15	11	17
2008 milpas	8	6	2	14	8	16

; and Supplementary Materials Tables S1 and S2). Of the 503 possible plazuelas identified using the TPI raster image, 111 aligned with previously documented plazuelas from the UAP survey, and 385 were newly identified possible plazuelas (Table 1). The TPI had a 35.2% confirmed positive rate through the identification of 111 of the 315 previously surveyed plazuelas. The plazuelas identified in the TPI were also compared to the 68 hilltops that have no known archaeological features, based on the UAP survey (Supplementary Materials Table S3). The TPI identified seven (10.3%) false positives. Because only 35% of known plazuelas were identified, it is possible that the 385 new possible plazuelas represent only a fraction of the actual plazuelas located within the UAP lidar zone. Overall, the output of the TPI proved to be useful for remotely identifying possible Classic Maya plazuelas.

3.2. Simple Local Relief Model (SLRM)

I used the open-source RVT [30,68,71] to create the SLRM raster model. Overall, the SLRM produced the best model of the three relief visualizations and OBIA. Classic Maya plazuelas are most visible on the SLRM (Figure 4f, Table 1 and Table 2). The SLRM allows the user to see the modified hilltops of possible plazuelas, as well as possible structures, based on small differences in local topography, due to human modifications to the landscape. These include transitions in topography where the sides of building platforms meet the plaza and the edge of the plaza at the crest of the hillslope.

In total, I identified 580 possible plazuelas across the UAP lidar dataset (Figure 6 and Table 1). Of the 580 possible plazuelas, 381 were also identified using TPI, and 199 were only identified on the SLRM (Table 2). Among the 580 possible plazuelas, 117 (37.1%) align with confirmed positive plazuela groups from the UAP survey (Figure 7). Thus, 457 new possible plazuelas were identified using the SLRM. The possible plazuelas were compared to the 68 hilltops known to contain no archaeological features resulting in a false positive rate of 8.8% (n=6) within the previously surveyed area (Figure 7). SLRM identified approximately 37% of known plazuelas, suggesting that hundreds more possible plazuelas exist in addition to the 457 SLRM possible plazuelas. In the upland landscape of the study region, SLRM was the most useful model for archaeological prospection.

3.3. SLRM and TPI Comparisons to Known Archaeological Features

The size of plazuelas and number of structures within each plazuela was compared between all surveyed plazuelas (n = 315) and confirmed plazuelas identified in remote archaeological prospection (n = 129). SLRM identified 117 previously surveyed plazuelas and TPI identified 111 previously surveyed plazuelas, 99 of which overlapped. TPI identified 12 surveyed plazuelas that were not identified by using SLRM, and SLRM identified 18 surveyed plazuelas that were not identified by using TPI. In total, 129 previously surveyed plazuelas were identified by using the two relief visualization techniques; of these 129 plazuelas, the plazuela area is known for 127, and the number of structures is known for 125 plazuelas.

3.3.1. Plazuela Size and Vegetation Cover

In general, archaeological prospection identifies larger plazuelas in terms of the size of plazuela and the size and number of structures in the plazuela group. However, this relationship is not consistent and is also impacted by vegetation. Previously surveyed plazuelas vary in size, from 43 to 11,760 m², although of the 315 plazuelas, 307 are less than 5000 m², and eight are larger than 5000 m² (Supplementary Materials Table S2). The largest plazuelas are the civic and ceremonial core areas of Uxbenká, which underwent massive anthropogenic landscape modifications during the first half of the Early Classic (250–600 CE) [52,76,77]. The average area of the previously surveyed plazuelas is 976 m² (

Table 3. Size of plazuelas and number of structures (based on pedestrian survey data) that were identified with SLRM and TPI. Larger plazuelas are more likely to be remotely identified. * Some plazuelas do not have areas calculated or number of structures recorded.

	Average Area (m2)	Sample Size *	Average Number of Surveyed Structures	Sample Size *
All previously surveyed plazuelas	976.3	315	2.8	305
SLRM only	943.3	17	3.4	17
SLRM and TPI	1613.6	98	3.5	97
TPI only	855.6	12	3.1	12
All SLRM	1514.5	115	3.5	114
All TPI	1530.9	110	3.5	109
Total all SLRM and TPI plazuelas	1452.3	127	3.5	125
All previously surveyed plazuelas (except plazuelas > 5000 sq. m)	804.1	307	-	-
SLRM and TPI (except plazuelas > 5000 sq. m)	1127.4	91	-	-
All SLRM (except plazuelas > 5000 sq. m)	1098.4	108	-	-
All TPI (except plazuelas > 5000 sq. m)	1095.7	103	-	-
Total all SLRM and TPI plazuelas except 7 largest (> 5000 sq. m)	1074.1	120	-	-

).

The remotely identified plazuelas vary in size from 223 to 11,760 m². Seven of the eight plazuelas greater than 5000 m² were remotely identified. However, no plazuelas under 223 m² were identified by using SLRM or TPI. These smaller plazuelas account for approximately 11% (n = 36) of the 315 surveyed plazuelas. The average size of plazuelas identified through remote sensing was higher when the disproportionately large plazuelas (greater than 5000 m²) were included (Table 3). These results varied from 1452 m² (total of all SLRM and TPI) to 1613 m² (plazuelas identified with both SLRM and TPI) (Table 3). Even when the plazuelas greater than 5000 m² are removed, the average size of remotely identified plazuelas remains larger than the average of the 307 plazuelas (804 m²). These results varied from 1074 m² (total of all SLRM and TPI) to 1127 m² (plazuelas identified with both SLRM and TPI) (Table 3). In general, larger plazuelas were identified by both SLRM and TPI, while smaller plazuelas varied, being identified by only SLRM, only TPI, or both techniques (Figure 8).

In addition to the size of the plazuela, vegetation cover impacts the ability to remotely detect plazuelas. Even small modified hilltops are visible if the hilltop was a milpa in 2011 when the lidar was acquired (Figure 8f). Medium plazuelas varied in their visibility, with densely vegetated plazuelas being more difficult to detect than plazuelas in 2011 milpas (Figure 8e,f). The largest plazuelas were visible even with dense forest cover (Figure 8d).

3.3.2. Number of Structures and Structure Size

Plazuelas with more structures and larger structures are more likely to be identified by using SLRM and TPI. Generally, larger plazuelas contain more structures and larger structures than smaller plazuelas, resulting in greater visibility on the SLRM and TPI of larger plazuelas and associated architecture. However, this relationship is not consistent (Figure 8b,e). Structure counts among surveyed plazuelas vary from one structure to 12 structures (Supplementary Materials Table S2), with an average of 2.8 structures per plazuela (Table 3). Nearly a third of surveyed plazuelas contain only one structure (

Table 4. Total structure counts of surveyed plazuelas compared to structure counts of plazuelas identified with SLRM and TPI. Plazuelas with more structures are more likely to be remotely identified.

Number of Structures	Surveyed Plazuelas Containing the Number Structures in the First Column	SLRM and TPI Plazuelas Containing the Number Structures in the First Column	% Identified
unknown	10	1	10%
1	100	24	24%
2	60	20	33%
3	53	23	43%
4	39	24	62%
5	26	14	54%
6	13	11	85%
7	7	5	71%
8	4	2	50%
9	1	1	100%
10	1	0	0%
12	1	1	100%
Total plazuelas (not including unknowns)	305	125	41%

). Across all surveyed plazuelas that were identified with SLRM and TPI, confirmed structure counts also varied from 1 to 12 with a 3.4 average number of structures. However, plazuelas with higher structure counts were remotely identified more often, with 24% of all plazuelas with one structure identified with SLRM and TPI, whereas any plazuela with more than four structures had a greater than 50% confirmed positive rate (Table 4).

These results parallel the findings from the previous study that examined a 4 km² subsample of the UAP lidar [9] but also indicate the advantages to using multiple relief visualization techniques. The previous 2015 study concluded that structures were difficult to detect on the bare earth hillshade and in the point cloud, and that hilltop modifications were most visible in the varied topography and vegetation [9]. Here, the SLRM and TPI outputs are useful for remotely identifying plazuelas; however, there is a bias for larger plazuelas and plazuelas with more structures.

3.4. Archaeological Prospection Compared to Land Use

Vegetation height and density varies based on land use practices. Previous studies observed that vegetation type and height across the landscape has a direct impact on the visibility of archaeological remains [6,9,35,41,62,78]. Modern Maya populations in Southern Belize practice slash-and-burn (swidden) agriculture resulting in a mosaicked patchwork of forest re-growth [9,79,80]. Previously, using a combination of a false color infrared (FCIR) GeoEye-1 multispectral satellite imagery to pinpoint color variations among patches of vegetation regrowth, lidar Digital Surface Model (DSM) which reflects vegetation foliage, and LAS point clouds, milpas from 2008 to 2011 in a 40 km² area within the UAP lidar dataset were mapped. These mapped milpas were used to test if vegetation growth (as a direct reflection of which year the land patch underwent milpa cultivation) influenced the visibility of archaeological remains using the SLRM and TPI model and comparing them to surveyed plazuelas.

milpas cleared in spring 2011 had an average of 6.3 ground returns per m². milpas cleared in 2010 to 2008 had decreasing ground returns due to the re-growth of tropical vegetation during fallow periods. In 2010 and 2009 milpas, ground returns were 3.9 and 3.3 returns per m². In 2008 milpas, the average ground returns of 1.6 per m² represented a nearly 75% reduction, as compared to 2011 milpas [9] (pp. 5–7). Recent milpas from 2009 to 2011 have higher average ground returns than milpas older than 2008 and the associated vegetation growth.

An example of the impact of milpas is visible at Settlement Group 83 from Uxbenká. Settlement Group 83 is divided into two plazuelas, A and B, which are similar in size: Plazuela A is 508 m², and Plazuela B is 442 m². While these plazuelas are of similar size, Plazuela A, which is in a 2011 milpa, is visible on both the SLRM and TPI, but Plazuela B, which is in a more vegetated patch of forest, is not visible on the SLRM or TPI. The small structures are not easily discernable in either Plazuela A or B (Figure 9).

Within the boundary of mapped milpas, 36.8% of confirmed plazuelas were identified using SLRM and TPI (

Table 5. Surveyed, SLRM, and TPI plazuelas identified within milpas from 2008 to 2011 (see Reference [9]). (* Note: One SLRM and TPI plazuela identified within the milpa boundary is a false positive).

	milpa Year				Total in the Mapped milpa Boundary
	2008	2009	2010	2011
Total plazuelas through SLRM and TPI	16	17	20	54	229 *
Possible plazuelas through SLRM and TPI	13	12	18	41	165
Confirmed plazuelas through SLRM and TPI	3	5	2	13	63
Confirmed plazuelas total	13	18	8	20	171
Percent of confirmed plazuelas identified with SLRM and TPI	23.1	27.8	25.0	65.0	36.8

). These results parallel the findings of the larger dataset, with a total of 129 confirmed positives among the 315 confirmed plazuelas (41%). Within the annual milpas, 65% of surveyed plazuelas were confirmed positives within 2011 milpas. The ability to remotely detect archaeological features, including modified hilltops, drastically decreases within a year of vegetation regrowth. Among 2008 and 2010 milpas, the confirmed positive rate was reduced to roughly 25% (Figure 10 and Table 5). Ultimately, the ability to remotely identify archaeological features varies based on the plazuela size, structure size, and vegetation cover, and while general trends are present (larger plazuelas with bigger structures are easier to detect), these relationships are inconsistent among medium- and smaller-size plazuelas, due to variations in vegetation cover.

4. Discussion

This study tested three lidar relief visualization techniques and OBIA on a slope model for the identification of archaeological features—plazuelas and house mound structures—in Southern Belize. Using a lidar-derived DEM, I used OBIA (on slope), SVF, TPI, and SLRM to assess remote archaeological prospection of hinterland Classic Maya settlements. Ultimately, TPI and SLRM were used to identify possible plazuelas. These new data were compared with pedestrian survey data for both the location of plazuela groups, the location of hilltops that do not contain archaeological features, and the impact of modern land use on archaeological prospection. These findings illuminate how lidar relief visualization techniques vary for archaeological prospection in the context of the southern foothills of the Maya Mountains in Southern Belize.

One of the largest challenges in Maya archaeology is the identification of hinterland settlements. While satellite imagery is a cost-effective way to remotely conduct archaeological prospection, the major limiting factor of freely available satellite imagery, including the newest versions, such as Landsat-8 and Sentinel radar imagery, is the relatively low resolution compared to the size of hinterland settlements. Additionally, because satellite imagery is a passive method of remote sensing (whereas lidar is an active remote sensing method), the vegetation in satellite imagery prevents and obscures visual access to the bare ground. Satellite imagery is useful for the detection of large archaeological sites and associated administrative architecture [5,81], and, while the detection of small hinterland (i.e., non-elite) settlements and household plazuelas has remained elusive [47], analyses of ecological zones were used to model potential hinterland populations [82]. The use of satellite imagery for archaeological prospection is complemented by the lidar revolution in archaeology since 2000.

The seminal lidar survey at Caracol in 2009 proved the potential of lidar-derived products for modeling ancient landscapes and detecting hinterland settlements in the Maya region [83]. Across the Maya Lowlands, archaeologists have relied on a variety of relief visualization techniques to test the applicability of lidar for the identification of archaeological features [1,4,8,9,33,35,41,43,59,84]. Combined, these previous studies suggest that no single technique is best for archaeological prospection due to the varied geography, including topography and modern human land use, and character of archaeological features such as the height/elevation, depth, size, and geometrical shape. In Southern Belize, SLRM and TPI were the most useful of the four methods evaluated for the identification of Classic Maya settlements.

Previous studies successfully applied OBIA and segmentation to lidar data [20]. For example, OBIA combined with Automated Feature Extraction (AFE) was used for the detection of archaeological sites in the Southeastern United States [21,25] and to identify more than 10,000 grave mounds in Tonga [44]. In Southern Belize, however, the homogenously hilly landscape resulted in a less-than-ideal OBIA raster for the detection of archaeological features. This is likely because the variations in topography are visually similar to the minor topographic changes from Classic Maya structures and plazuelas. OBIA on a lidar slope model would likely work well in a location that is relatively flat, as the topographic images would “pop out” or be identified with segmentation methods.

Similar to the studies with OBIA, other studies have successfully used SVF in archaeological prospection. Scholars studying prehistoric hillforts in Slovakia [29] and Slovenia [30] used SVF for the identification of archaeological features. SVF was successfully applied to the detection of ancient Maya plazuelas and structures at El Pilar [17]—which is a protected forest reserve—as well as Chichén Itzá [35]. SVF proved less successful both on the UAP data and in the Northern Maya Lowlands at Ucanha [6]. The quality of the SVF output for identifying archaeological sites is not likely caused by the direction input features (see Reference [71]) but is more likely due to the variations in vegetation height.

In this study, the two most successful relief visualization techniques were SLRM and TPI. Combined, 563 new possible plazuelas were identified in this study. Other archaeological studies in the Maya region, including the Belize River Valley [41] and Caracol [40], used TPI to identify archaeological features including terraces and house mounds. In the Belize River Valley, TPI was used to identify ancient Maya house mounds among a modern semi-urban community and to test the impact of land use (developed vs. undeveloped) and vegetation class (pasture, orchard, or forest) on the visibility of archaeological remains using a TPI landscape model [41]. These techniques have also been used outside of the Maya region, in places such as Western Romania, for the detection of hillforts [38].

Here, I identified a total of 702 possible plazuelas using both TPI and SLRM. Out of the 315 surveyed plazuelas, only 41% (n=129) were identified by using SLRM and TPI (Table 2). The confirmed positive rate (41%) of identified settlements on the UAP lidar dataset is similar to other confirmed positive rates in the Maya region. In the Yaxuná-Popola-Tzacauil area, 32% of mapped residential structures were identified in lidar relief visualizations [35]. In Western Belize, 27% (4 of 15) mapped residential structures were identified on the lidar hillshade and TPI raster [42]. In the Northern Yucatan, 47% of mapped buildings were visible on the lidar between Ucí and Cansahcab, and at Ucanha, 48% of architectural features were visible [7].

In total, 563 new possible plazuelas were identified (702 total plazuelas, minus 129 surveyed plazuelas and 10 false positives). If the 563 new possible plazuelas represent a similar proportion of remotely identified plazuelas (approximately 40%), it is possible that more than 1370 plazuela groups are present on the landscape within the UAP lidar zone. Factoring in a 14.7% rate for false positives, the number of new possible plazuelas is still over 1150. The ability to remotely detect large plazuelas informs the extent of the settlement system, but the inability to detect more than half of the surveyed plazuelas groups, and especially smaller plazuelas groups often associated with hinterland households [49], highlights the need for ground-truthing in conjunction with remotely identified archaeological features [17,35].

Previous settlement pattern studies suggest diversity in the size, clustering, and distribution of Classic Maya plazuelas in Southern Belize [50,85]. The diversity and variability in the size of documented plazuelas is likely representative of the variability in the size of newly identified plazuelas. Furthermore, in the Maya region, the size of households—in this study, they are referred to as plazuelas—is often directly linked to land tenure and the intergenerational transmission of wealth, where the oldest households develop into the largest and wealthiest households [49,86]. Others have used lidar to assess changing settlement patterns over time, based on temporal variations in architecture [4,43,84]. Among the newly identified plazuelas, larger plazuelas were distributed across the landscape, possibly indicating long-standing households in the hinterlands between Uxbenká and Ix Kuku’il and the Maya center Lubaantun to the east, although ground-truthing and test-unit excavations to gain chronologic information are needed to confirm this hypothesis. Little pedestrian survey has occurred between the regional Maya centers of Southern Belize, with the exception of a transect survey between Nim Li Punit and Xnaheb [57] and extensive survey between Ix Kuku’il and Uxbenká. The survey between Nim Li Punit and Xnaheb suggests a decrease in settlement density acting as a boundary between the centers [87]; similar trends have been noted at Uxbenká and Ix Kuku’il based on decreased settlement density in a river valley between the centers [48]. However, the results from this study suggest a more densely populated landscape between Lubaantun and Uxbenká, highlighting the importance of ground-truthing and diversity in settlement patterns and household clustering within the Southern Belize region [50].

Areas where archaeologists work are constantly changing as living landscapes. Within the UAP lidar zone, there are expanding modern villages; shifting agricultural fields where maize, beans, pumpkins, and rice are grown; tree orchards of cacao and palm groves; and, more recently, cow pastures [80,88,89]. These constantly changing diverse settings must be accounted for during the use of remote-sensing imagery to assess the landscape for archaeological features. For example, modern house building often requires leveling of the land, resulting in rectangular cuts in the landscape that are visible in the lidar data and resulting outputs (Figure 11). Likewise, modern footpaths, tracks (unpaved roads), and football (soccer) fields are visible on the lidar DSM, DTM, and TPI and SLRM raster files (Figure 11). The untrained eye could easily mistake these modern landscape modifications for ancient modifications associated with archaeological features. Furthermore, modern anthropogenic landscape modifications must be considered if conducting AFE [90].

Vegetation type, height, and density impact the lidar returns. In the Yucatan [6], Pasion region [8,43], and in Southern Belize [9,18], medium-growth vegetation results in decreased visibility of archaeological remains. The findings presented here and elsewhere [6,43] indicate that recently burned milpas and low pastures positively impact the ability to detect small structures. Here, the number of ground returns is higher for 2011 milpas than for older milpas or areas of forest regrowth (Figure 12 and Figure 13).

In addition to the height of vegetation, the type and density of vegetation impacts archaeological prospection [10,62], especially in the identification of small architectural features such as house mounds (see Figure 9). Within medium to high canopies, the variation in underbrush height and density impact the ability of lidar laser shots to hit the ground [6], resulting in patches of land containing little-to-no ground points [8]. These patches with few ground points directly impact the ability to remotely detect small archaeological features. Areas with recent milpas and low grass of cow pastures contain higher ground points than areas of forest regrowth or even orchards, where the dense, low foliage of the trees prevent the laser from passing to the ground (Figure 13b).

5. Conclusions

The use of remote sensing and lidar technology has revolutionized the effectiveness and efficiency of archaeological prospection and survey [11]. Every relief visualization technique has benefits and disadvantages based on topography, vegetation, and archaeology type, including feature size and shape [64]. In this case study from Southern Belize, lidar-derived relief visualization methods SLRM and TPI were selected to visually identify Classic Maya plazuelas and structures. In total, 580 possible plazuelas were identified by using SLRM, and 503 were identified by using TPI, of which 381 overlap. Combined, a total of 702 possible plazuelas were remotely identified, 129 of which align with previously surveyed plazuelas, and 10 were known false positives, for a total of 563 newly identified possible plazuelas. These data help shed light on the density and dispersal of Maya settlements in areas that lack pedestrian survey.

While the results of these methods underestimate the number of plazuelas, in conjunction with previous pedestrian survey data, it is possible to estimate the total number of Classic Maya plazuelas within the lidar area. In total, 563 new possible plazuelas were identified by using SLRM and TPI, although many more likely exist. The 41% confirmed positive rate of previously surveyed plazuelas indicates that upwards of 1370 previously unrecorded Classic Maya plazuelas could be present on the landscape. Considering the 14.7% false-positive rate of SLRM and TPI, combined there could still be more than 1150 possible plazuelas on the landscape that are yet to be documented. These results emphasize the need for ground-truthing, in order to identify false positives, record the number of the structures which are difficult to detect on the relief visualization, and survey areas that were not identified through the use of SLRM and TPI, as they may contain ephemeral archaeological features.

Manual feature extraction is effective and allows the user to apply checks and balances to their processes [19]. Here, I used previously surveyed data to assess the accuracy of remote archaeological prospection to the location and size of known plazuelas. Future work will expand on this study by using AFE to test the utility of machine intelligence for archaeological prospection, in comparison to manual identification of archaeological features. AFE is increasingly popular for archaeological prospection [91] and has been successful in a variety of archaeological contexts using satellite imagery [5,92,93]. The application of AFE to lidar data has taken off in recent years (see References [21,64,94]), but it is yet to be extensively applied to settlement studies for the identification of small archaeological features within the Maya Lowlands (see Reference [90] for a recent application of AFE on Airborne Laser Scanning data in the Maya region). Applications of AFE to settlement studies in the Maya region has the potential to further advance the geospatial revolution in Mesoamerican archaeology.

Relief visualization methods are useful for archaeological prospection, but results can vary based on local vegetation, landscapes, and type of archaeological features. Modern human behaviors such as expanding villages, agroforestry (orchards), cow pastures, and shifting agriculture cycles impact lidar data and the ability to remotely detect archaeological features in the different ways. Remote sensing and lidar relief visualization techniques such as SLRM and TPI can increase our understanding of Classic Maya settlement systems, forcing us to re-evaluate the extent of settlement systems, population estimates, and the relationships of people in the past.