Classic Maya civilization in detail Lidar (a type of airborne laser scanning) provides a powerful technique for three-dimensional mapping of topographic features. It is proving to be a valuable tool in archaeology, particularly where the remains of structures may be hidden beneath forest canopies. Canuto et al. present lidar data covering more than 2000 square kilometers of lowland Guatemala, which encompasses ancient settlements of the Classic Maya civilization (see the Perspective by Ford and Horn). The data yielded population estimates, measures of agricultural intensification, and evidence of investment in landscape-transforming infrastructure. The findings indicate that this Lowland Maya society was a regionally interconnected network of densely populated and defended cities, which were sustained by an array of agricultural practices that optimized land productivity and the interactions between rural and urban communities. Science, this issue p. eaau0137; see also p. 1313

Structured Abstract INTRODUCTION Lowland Maya civilization flourished from 1000 BCE to 1500 CE in and around the Yucatan Peninsula. Known for its sophistication in writing, art, architecture, astronomy, and mathematics, this civilization is still obscured by inaccessible forest, and many questions remain about its makeup. In 2016, the Pacunam Lidar Initiative (PLI) undertook the largest lidar survey to date of the Maya region, mapping 2144 km2 of the Maya Biosphere Reserve in Guatemala. The PLI data have made it possible to characterize ancient settlement and infrastructure over an extensive, varied, and representative swath of the central Maya Lowlands. RATIONALE Scholars first applied modern lidar technology to the lowland Maya area in 2009, focusing analysis on the immediate surroundings of individual sites. The PLI covers twice the area of any previous survey and involves a consortium of scholars conducting collaborative and complementary analyses of the entire survey region. This cooperation among scholars has provided a unique regional perspective revealing substantial ancient population as well as complex previously unrecognized landscape modifications at a grand scale throughout the central lowlands in the Yucatan peninsula. RESULTS Analysis identified 61,480 ancient structures in the survey region, resulting in a density of 29 structures/km2. Controlling for a number of complex variables, we estimate an average density of ~80 to 120 persons/km2 at the height of the Late Classic period (650 to 800 CE). Extrapolation of this settlement density to the entire 95,000 km2 of the central lowlands produces a population range of 7 million to 11 million. Settlement distribution is not homogeneous, however; we found evidence of (i) rural areas with low overall density, (ii) periurban zones with small urban centers and dispersed populations, and (iii) urban zones where a single, large city integrated a wider population. The PLI survey revealed a landscape heavily modified for intensive agriculture, necessary to sustain populations on this scale. Lidar shows field systems in the low-lying wetlands and terraces in the upland areas. The scale of wetland systems and their association with dense populations suggest centralized planning, whereas upland terraces cluster around residences, implying local management. Analysis identified 362 km2 of deliberately modified agricultural terrain and another 952 km2 of unmodified uplands for potential swidden use. Approximately 106 km of causeways within and between sites constitute evidence of inter- and intracommunity connectivity. In contrast, sizable defensive features point to societal disconnection and large-scale conflict. CONCLUSION The 2144 km2 of lidar data acquired by the PLI alter interpretations of the ancient Maya at a regional scale. An ancient population in the millions was unevenly distributed across the central lowlands, with varying degrees of urbanization. Agricultural systems found in lidar indicate how these populations were supported, although an irregular distribution suggests the existence of a regional agricultural economy of great complexity. Substantial infrastructural investment in integrative features (causeways) and conflictive features (defensive systems) highlights the interconnectivity of the ancient lowland Maya landscape. These perspectives on the ancient Maya generate new questions, refine targets for fieldwork, elicit regional study across continuous landscapes, and advance Maya archaeology into a bold era of research and exploration. Representation of the archaeological site of Naachtun, Petén, at twilight. Each ancient structure is marked by a yellow dot. CREDIT: L. AULD-THOMAS AND M. A. CANUTO

Abstract Lowland Maya civilization flourished in the tropical region of the Yucatan peninsula and environs for more than 2500 years (~1000 BCE to 1500 CE). Known for its sophistication in writing, art, architecture, astronomy, and mathematics, Maya civilization still poses questions about the nature of its cities and surrounding populations because of its location in an inaccessible forest. In 2016, an aerial lidar survey across 2144 square kilometers of northern Guatemala mapped natural terrain and archaeological features over several distinct areas. We present results from these data, revealing interconnected urban settlement and landscapes with extensive infrastructural development. Studied through a joint international effort of interdisciplinary teams sharing protocols, this lidar survey compels a reevaluation of Maya demography, agriculture, and political economy and suggests future avenues of field research.

Lowland Maya civilization flourished for nearly 2500 years (1000 BCE to 1500 CE) in southern Mexico, Guatemala, and Belize (Fig. 1), the central part of which consists of approximately 95,000 km2 of rolling karst topography interspersed with wetlands (1). Today, this area is largely covered by tropical forest that has severely limited the scale of on-the-ground archaeological research, hindering regional assessments about ancient urbanism, population size, resource management, and sociopolitical complexity.

Fig. 1 Distribution of 12 PLI survey blocks. Location of 12 PLI survey blocks in relation to the central Maya Lowlands area as well as to the sites mentioned in text; for those survey blocks without a named site, the name of the block is provided (i.e., Env 1, Env 2, and Yala).

In 2016, the Pacunam Lidar Initiative (PLI) undertook the largest lidar survey to date of the lowland Maya region, mapping a total of 2144 km2 of the Maya Biosphere Reserve (MBR) in Petén, Guatemala. This survey mapped 10 noncontiguous areas (polygons) ranging in size from 91 to 454 km2, making it possible to characterize ancient settlement, agriculture, and physical infrastructure over an extensive and varied swath of the central Maya Lowlands by using a standardized and consistent dataset. Previous lidar surveys had taken place in Belize, Guatemala, and Mexico (2–8), providing important data on ancient lowland Maya civilization. However, the PLI sample, undertaken at far larger scale, captured greater variety in topography and ancient settlement for the central Maya Lowlands, the demographic and political heartland (9) of the Classic Maya culture. These data provide an opportunity to address current debates in Maya archaeology and to deepen studies of complex, tropical civilizations in the ancient world.

Some scholars suggest that the Maya Lowlands contained small city-state centers ruled by warring elites (10–13), in settlements supported by a relatively sparse rural population practicing swidden farming (14), with only limited input from intensive agriculture (15–21). In contrast, other views point to a regional network of densely populated cities with complex integrative mechanisms (22–25) that depended on heavy labor investments inside and outside urban cores (15, 26, 27). Even though the latter view has been ascendant in recent years, the absence of regional data has left the debate unresolved. The PLI data cover a sufficiently large area to provide robust support for the latter model and offer further insights about human-environment interactions in the region, including unexpectedly extensive fortifications and road networks.

Data collection, processing, and field validation Lidar is an active remote sensing technology that uses laser pulses to map landcover and ground surface in three-dimensional space. The PLI effort collected lidar data with a Teledyne Optech Titan MultiWave multichannel, multispectral, narrow–pulse width lidar system, the most advanced aerial lidar system (28) deployed in the region to date (Table 1). The improvements include scanning of the terrain from six different view angles (versus the more usual two views) and higher sensor-range resolution, enabling more precise mapping of topographic signatures under short vegetation. The data collection was designed with a nominal laser density of 15 pulses/m2 from a flying altitude of 650 m above ground level. The fidelity of this sensor reveals details of local topography, ancient architectural and agricultural features, and damage from looting (Fig. 2). Table 1 Lidar ground point averages. Important archaeological sites and lidar ground point densities within each PLI polygon. View this table: Fig. 2 Close-up of Xmakabatun. Image of Xmakabatun that demonstrates the high-fidelity detail of the Teledyne Opech Titan sensor. (A) Looting depressions highlighted by an openness visualization. (B) Looting features as drawn in the field (shown in red) as well as evidence of monumental architecture, causeways, residential structures, ditches, and terraces. The data were collected in 10 polygons over 12 designated survey blocks (Table 1). Each survey block was named for the largest known archaeological site within its boundaries; two survey blocks in archaeologically unknown areas were designated Env 1 and Env 2. The current dataset arose from roughly 33.5 billion laser pulses, producing about 60 billion returns, of which 5.2 billion were classified as ground returns. The average ground densities for different polygons range from 1.1 to 5.3 returns/m2, the result of varying vegetation densities and canopy coverage across the survey area. Of the total mapped area, 23% has an average ground return density of 1.1 returns/m2, 60% has a ground return density range of 1.5 to 2.7 returns/m2, and 17% has an average ground return density of more than 5.2 returns/m2. These data were interpolated to create bare-earth terrain models with a spatial resolution of 1 m per grid cell. Identification of archaeological features within the lidar dataset was iterative (29). The PLI study region has been subject to thousands of person-days of pedestrian survey over the course of eight decades of field research. Consequently, there was an extensive “ground-truthed” dataset available for most survey blocks (Fig. 3). Analysts from each participating team made preliminary identifications using existing site maps as training samples. Features of archaeological interest were identified through visual inspection of lidar-derived terrain models. Multiple visualization rasters and manual identification were used to highlight distinct kinds of topographic features (Fig. 4); Table 2 specifies the range of greatest utility in such visualizations. In addition to previously published methods (30–32), we developed a new technique, “prismatic openness” (33), that facilitated detection of small features such as simple house platforms and agricultural terraces. Fig. 3 Mapped areas within PLI zone. Map of survey blocks showing (i) areas covered by pedestrian surveys before acquisition of lidar data and (ii) areas where ground verification of the lidar data has been completed to date. Fig. 4 Visualizations used for lidar analysis. (A) Red Relief Image Map (RRIM); (B) Sky-View Factor (SVF); (C) Simple Local Relief Model (SLRM); (D) Prismatic Openness. Table 2 Lidar visualizations. Raster visualization combinations created and used for feature identifications. View this table: After preliminary identifications backed by preexisting map and excavation data, ground-truthing and test excavations occurred in 2017 within the Holmul, Corona, Perú, Naachtun, Tintal, Uaxactun, and Zotz survey blocks. This allowed each project to (i) calibrate their feature identifications to ground conditions and (ii) identify patterns in false-positive and false-negative identifications. With this knowledge in hand, a second wave of feature identification took place from August to December 2017. Overall, approximately 165 km2 within the study area have either been used as training samples, subsequently “ground-truthed,” or both (Table 3); that is, 7.7% of the entire PLI survey region has been verified through ground survey, although mapped areas were not uniformly distributed (Table 3). Moreover, all field teams within the PLI study region have ground-validated features belonging to every class reported on here: buildings of various types, defensive features, upland and wetland agricultural features, causeways, canals, and reservoirs. Table 3 Extent mapping and ground verification. Areas in the PLI survey region covered by both pre-lidar pedestrian surveys and post-lidar ground validation efforts (these sometimes overlap). View this table: PLI’s preliminary field validation efforts suggest that the feature identifications are accurate if slightly conservative. Most projects report a net increase in field-verified structures of ~15% over those identified through visual analysis of the lidar data. Furthermore, certain features—in particular, small channels or berms as well as large causeways—are clearly visible in the lidar although almost impossible to recognize on the ground without excavation, meaning that in-field visual inspection alone is not always capable of assessing the accuracy of lidar feature identifications. Nonetheless, because ground-truthing efforts are ongoing, and because error rates (false positives and negatives) are likely to vary for different classes of features, all analyses presented here use the 2017 PLI feature dataset with no further adjustments.

Population estimates Using the above methodology, 61,480 structures have been identified in the PLI survey region, resulting in an aggregate settlement density of 29 structures/km2 (Table 4). Because this structure count derives from a much-expanded and regionally representative dataset, it can serve as the basis for a uniquely robust regional population estimate. Even so, producing population estimates from archaeological data requires some assumptions of varying detail and scope from ethnographic analogy and archaeological data; as a consequence, the exercise comes with some degree of uncertainty that only future field research can reduce. Table 4 Total structures and settlement density estimates by PLI survey blocks. View this table: To estimate Late Classic population from structure counts, well-established methodologies for the Maya Lowlands were adjusted by controlling for invisible structures, contemporaneity of use, nonresidential function, non–Late Classic occupation, and numbers of individuals per building (Table 5) (34). The values of these adjustments vary widely across the literature, given the regional variability of the archaeological record. The geographic scope of the PLI study area, which includes diverse landscapes with distinct settlement histories, necessitates a broader view beyond any single set of local adjustments. As such, values for each adjustment were collated from published scholarship and categorized into five sets according to their overall impact (minimum, low, middle, high, and maximum). All adjustments were then multiplied to derive a composite population index range (PopIndex), which, when multiplied by the total number of structures, produced an overall population range (29). Table 5 PopIndex calculations. Set of adjustments used to calculate a “population index” (PopIndex), which produces an estimated population when multiplied by structure number. The range in the value of each adjustment indicates the variation within the published record. View this table: Population estimates were calculated to the nearest half million to avoid false precision and are presented as a range rather than a single value for the central Maya Lowlands. Given the nature of the available excavation data, this study does not propose a diachronic demographic curve. It only estimates total population during the Late Classic population maximum, 650 to 800 CE, because all PLI survey blocks subjected to archaeological investigation have shown evidence for substantial Late Classic populations. Following this procedure and guidelines, we estimate a Late Classic population range of 150,000 to 240,000 for the entire PLI survey region. This number amounts to an average density of ~80 to 120 persons/km2. This density value accords with the ~100 persons/km2 suggested by previous nonlidar research (9, 34–37). Inomata et al.’s recently published lidar-derived settlement data from the Ceibal region (with 31 observed structures/km2 across a sample of 470 km2) (38) are consistent with the settlement data presented here. There is substantial variability in settlement density across the PLI survey blocks, which reflects a pattern characteristic of the entire central Maya Lowlands. That is, some areas will be densely settled (comparable to the Tikal, Xultun, or Naachtun survey blocks) and others will be sparsely occupied (comparable to the Corona, Yala, and Env 1 and 2 survey blocks). However, because average settlement density over large areas is likely to be similar, the PLI dataset can be extrapolated with confidence over a 95,000-km2 area known as the “central Maya Lowlands” (Fig. 1). This area excludes the northern Yucatan, the Gulf Coast plains, and coastal Belize and Quintana Roo because they differ physiogeographically from the PLI survey region. The Late Classic population for the central Maya Lowlands was calculated by multiplying (i) the mean structure density (structures per km2) of the sample, (ii) the total area of the central Maya Lowlands (95,000 km2), and (iii) the PopIndex range (see below): (1)This formula suggests that 61,480 identified structures within 2144 km2 of PLI’s surveyed area can be extrapolated to 2,724,396 structures over an area of 95,000 km2. Once multiplied by the PopIndex range, these numbers yield a population range of 7 million to 11 million for the central Maya Lowlands during the Late Classic period. Notably, this range aligns with the upper end of prior estimates for the same region, derived from pedestrian surveys (36, 39–41). The incorporation of Inomata et al.’s published data from Ceibal (38) to both PLI’s structure and coverage area totals results in the same extrapolated population levels for the central Maya Lowlands—a good indication that the PLI dataset is a representative sample for the region and that the extrapolations are sound. This estimate is moderate, even conservative, primarily for two reasons: (i) The traditional method of calculating regional population—multiplying area (95,000 km2) by a constant population density (80 to 120 persons/km2)—has been found to underestimate the resident populations of larger centers (Tikal, Naachtun) that become denser with increasing scale (42). Accordingly, the densest urban areas likely would have contained more people than these current approximations. (ii) Although comparison of lidar-derived structure data and pedestrian survey data shows strong agreement within the PLI survey area, vegetation and taphonomic factors reduce the visibility of small structures in lidar-based visualizations. Ground validation of the PLI lidar-derived dataset has resulted in a ~15% net increase in archaeological structures. Because both of these factors would result in larger populations, only further fieldwork can render these estimates more accurate. Regardless of these caveats, the PLI data suggest overall regional populations of such scale that some degree of agricultural intensification (43–45) would have been required to sustain them. Approximately half of the central Lowlands are seasonal wetlands known as bajos (15). Because permanent settlements tended to avoid these flood-prone and poorly drained areas, they remained largely uninhabited and could then be available, after added investment, to intensive agriculture (46).

Agricultural intensification Early applications of airborne synthetic aperture radar (SAR) suggested the existence of wetland field systems throughout the entire central Maya Lowlands (47–49). However, most of the patterns identified as evidence of intensive agricultural features were later dismissed as “noise” in the imagery or artifacts of the local geology (15, 21, 50–52). Furthermore, given that most of the wetlands in the interior of the peninsula contain clay vertisols, their overall suitability for ancient agriculture was challenged (17, 50, 52–57). Although aerial photography and fieldwork succeeded in identifying both upland and wetland intensive agricultural features in various zones within the lowlands (58–60), the extent of known agricultural improvements remained limited to a few specific areas (18–20, 55, 61–67). Even though the idea that lowland Maya society only relied on swidden was discarded in the 1980s in favor of a “new orthodoxy” claiming some reliance on a mosaic of more intensive agricultural strategies, the dominant view among specialists was that wetland agriculture was generally limited to areas peripheral to the central Maya Lowlands and that the elevated interior would have relied primarily on swidden agriculture (15, 18, 19, 27, 50, 52, 54, 61–63, 68, 69). The PLI survey revealed a landscape heavily modified for intensive agriculture: 11 of the 12 survey blocks include agricultural features of some kind. In wetland areas, the survey primarily identified drainage channels. These were designed either to draw water away from natural streams toward infrequently flooded areas or to drain those same areas during major floods. Their dimensions vary from 1 to 2 m in width and 20 to 50 cm in depth. Major channels can measure more than 1 km in length (one extreme case at Tintal is 2.5 km long). Large and small channels intersected at regular intervals, forming nested grids within “channelized” fields (Fig. 5). Overall, these data suggest that wetlands throughout the central Maya Lowlands were modified for agricultural exploitation (16, 43, 61, 62, 70), rather than being used exclusively as sources of water and transportable soils (16, 71). The largest total area of wetland agriculture is found in the Naachtun region, with 31.5 km2 of interspersed channelized zones, followed by Holmul and Tikal. The largest single zones of contiguous wetland fields occur in the Holmul survey block, ranging from 2.3 to 7 km2, followed by Tikal and Naachtun. These wetlands may have been selected because of local edaphic conditions, yet they also lie near the most populated areas. Fig. 5 Wetland and upland agricultural features. Examples of intensive wetland fields featuring networks of canals (above) and upland field zones featuring low stone wall features and terraces. On upland terrain, numerous linear stone features such as terraces and low field walls were identified. These features measure 1 to 2 m in width and are no higher than 50 cm. Some are strictly linear (i.e., contour terraces and check dams); others create enclosures (Fig. 5). Several features form gridded fields with or without enclosure walls, and, in some examples, surround residential structures similar to the albarradas (circular walls) of the northern Maya Lowlands (72). Because of their multiple functions, stone features vary widely in length, from only a few meters to >1 km. By identifying these features, we defined “zones of intensive cultivation” (29) where agricultural improvements were either the only feature in the landscape (e.g., wetland fields) or were interspersed with settlement. The latter case encompasses production systems generally described as infields, orchards, and houselot gardens. A total of 67 km2 (6659 ha) of wetland fields and 295 km2 (29,517 ha) of upland fields were classified as zones of intensive cultivation, accounting for ~17% (362 km2) of the entire PLI survey area. Investments in upland agriculture are most common in the central and eastern sectors of our sample. The largest area totals occur, by descending order, in the Tikal (64 km2), Uaxactun (57 km2), Naachtun (47 km2), and Xultun (45 km2) survey blocks. Estimates for the overall productive capacity of this heavily modified agricultural landscape were calculated to determine whether it could have supported the Late Classic populations reported above. Little ethnographic information on the productivity of intensive agriculture in the Maya Lowlands is available, so estimations of overall productive capacity must rely on the productivity of swidden agriculture as the sole point of reference. Because some of the improved agricultural land was unusable (e.g., bajo) or prone to rapid degradation (e.g., erosible slopes), many improvements would have succeeded in extending productivity to otherwise unproductive areas (73). However, because some upland features were indeed placed on optimal land (i.e., gentle slopes) and the moisture-retaining properties of terraces improve agricultural productivity (37, 74), it is likely that some upland features did improve overall yields beyond what traditional swidden agriculture could produce. Nonetheless, these calculations assume that all zones of intensive cultivation resulted in the same productive capacity as swidden. Ethnographic data on the productivity of traditional lowland Maya swidden agriculture were collated and standardized to estimate the number of hectares necessary to support one person. Data on swidden productivity were drawn from Cowgill (75–77), Griffin (78), Schwartz (79, 80), Nations and Nigh (81), and Ford and Nigh (37). Standardization entailed, among other adjustments, normalizing the annual productivity of different fallow regimes to multiyear averages. These calculations resulted in a mean productivity value of 0.48 ha per person [the mean of the range (0.34 to 0.62 ha per person) resulting from using different adjustments from the above-cited sources] (table S3). This value was then applied to all the land in the PLI sample (1314 km2) available for agricultural production (zones of intensive cultivation and nonwetland rural areas) to derive a maximum supportable population. Finally, by the same logic applied to population estimates, this number was adjusted to account for the possibility that these features of agricultural intensification were not all coeval. These calculations (Table 6) suggest that the total amount of land available for agricultural production in the PLI sample would have been more than sufficient to sustain the population levels reported above (29). Table 6 Maximum sustainable population of each survey block. Calculations designed to determine the maximum sustainable population within each survey block based on (i) its total measured agricultural area (values adjusted for Late Classic contemporaneity using AgroIndex) and (ii) a minimum swidden area of 0.48 ha/person (mean value of range reported in table S3). View this table: Consequently, the PLI data indicate some combination of the following possibilities: (i) There was capacity for surplus food production; (ii) substantial portions of agricultural land were dedicated to the production of nonfoodstuffs (such as cotton or cacao); (iii) large areas of the uplands remained covered in forest (37, 82, 83); or (iv) populations were higher than our current estimates. Agricultural intensification could also have been achieved using strategies such as fallow-period shortening, multicropping, and manipulation of crop sociology; these would not have occasioned changes to the landscape visible in lidar. In addition to corn, a variety of other edible plants and fruit trees are known to have been grown in and around improved upland fields (14). Therefore, this study’s reliance on both physical infrastructure (as an index for intensification) and ethnographic data of modern swidden productivity likely underestimates the total productive capacity of ancient lowland Maya agriculture. This dataset also demonstrates great variability in agricultural investment across the region (Fig. 6). The Corona survey block (432 km2) has <1% of its available land improved for cultivation, whereas the Tikal block (147 km2) has as much as 70%. Indeed, the lowest densities of agricultural features occur in the less populated western blocks of the sample, and the highest densities are located in the central and eastern sectors where settlement is densest. Across the PLI sample, areas defined here as periurban and urban (see below) contained the greatest relative density of agricultural features. Urban cores contained no agricultural features, whereas rural areas contained relatively fewer improved fields (Table 7 and Fig. 6). The scale of wetland field systems and their close association with densely populated areas suggest some degree of centralized involvement. Conversely, agricultural modifications of the uplands cluster around residential groups, which suggests that they were likely managed by household or small corporate groups. Upland improvements and terracing are common, likely making a substantial contribution to overall regional production (26). Fig. 6 Agricultural features relative to population densities. Map of survey blocks showing deviations from expected density of upland field areas relative to settlement density class area. Expected (i.e., “no difference”) value is 1.0. Table 7 Correlation of upland agricultural features with settlement density. Correlation matrix of intensive agriculture upland field zones and settlement density zones in m2 (upland zones strongly correlate with periurban and urban zones; N = 12). View this table: The co-occurrence of dense settlement with agricultural improvements suggests that investments to maximize agricultural production were directed to areas where population sprawl limited the possibility for extensive farming. Despite these efforts, however, the estimated yield of the agricultural lands within the most heavily modified survey blocks—Tikal, Naachtun, Xultun, and Perú—would not have been sufficient to sustain their populations (Table 6). Even if the zones of intensive cultivation within these survey blocks had been twice as productive as swidden fields, none of these areas would have harvested enough to achieve self-sufficiency. That is, people living in these blocks were partly dependent on foodstuffs imported from surrounding areas. Such data point to a lowland Maya landscape that was a mosaic of interdependent densely populated political centers and extensive periurban hinterlands engaged in a complex set of interactions.

Urbanization Settlement studies at lowland Maya sites have often used categories such as “site core,” “periphery,” and “rural” to delimit zones of different settlement density (84–86). These are useful categories at local scales. However, the scope of the PLI dataset allowed for a comparison of settlement density at a regional level. A shared settlement density scale was developed by applying figures compiled by Rice and Culbert (35) (Table 8). Their data, based on ground survey of sites with variably defined limits, provided the justification for four classes; a fifth—“vacant”—was added for a total of five classes. In summary, five structure density categories—vacant, rural, periurban, urban, and urban core—were defined for the entire dataset (29) to discern differences not only in the size and density of cities but also in the relationship between cities and hinterlands (86–89). Table 8 Late Classic structure density classes, modified from Rice and Culbert (35). View this table: The “urban core” class was defined on the basis of structure densities at the heart of the largest Maya centers, such as Tikal (~700 structures/km2). “Urban” describes the structure density at the heart of many smaller cities within the PLI study area (e.g., Uaxactun, Zotz) as well as the densely settled surroundings of major cities’ urban cores (e.g., Xultun, Naachtun, and Tikal). “Periurban” zones (also known in the literature as “urban fringes”) combine characteristics of both urban and rural areas (90). The parameters of this category encompass those areas designated by prior surveys as “peripheral” to the largest Maya sites as well as sites often characterized as “minor centers” whose economic and political fortunes were closely articulated with larger cities nearby (91–93). Our “rural” class follows existing descriptions of settlement density in the hinterlands of smaller cities and minor centers. Finally, we classified as “vacant” all zones with fewer than 10 structures/km2. At the regional level, sprawling urban and periurban settlement zones covered large areas in the east, whereas the west remained mostly rural (Fig. 7). Urban core densities ran as high as 1000 to 2000 persons/km2, whereas rural areas were as low as 50 persons/km2 over continuous areas. However, the PLI sample also revealed striking regional differences in urbanization patterns (Table 9). Fig. 7 Settlement density patterns. Density values are expressed as structures/km2. Table 9 Relative frequency of structures in each density class by PLI survey block. View this table: Cluster analysis of the settlement patterns within the PLI sample revealed several intraregional patterns. Because the PLI survey blocks varied in shape and in total area, standardized 80-km2 sample areas were designed for cluster analysis. These areas were delineated to maintain the primary urban center in as central a location as possible within the sample area. Because of the variable shape and dimensions of each PLI survey block, sample areas could not be drawn with identical proportions; nonetheless, areas were 80 km2 in all cases, with proportions as standardized as possible. The Corona and Holmul survey blocks were large enough to delineate two distinct sample areas for each survey block. These two additional sample areas were named Achiotal and Xmakabatun (for the Corona and Holmul survey blocks, respectively) after the name of the primary center within each area. Cluster analysis of these 14 areas (Fig. 8) (29) indicates that settlement patterns at the regional level fall into three classes: (i) rural areas with low overall density and limited political integration (Env 1, Env 2, Yala, Corona, Achiotal); (ii) periurban areas containing small urban centers with large periurban and rural populations (Uaxactun, Holmul, Xmakabatun, and Zotz); and (iii) urban areas with periurban and rural populations surrounding and likely integrated with a single large urban center. Within this third class are two subclasses: moderately urbanized areas with more substantial periurban and rural populations (Tintal and Xultun), and highly urbanized areas in which most of the population lived in urban core, urban, and periurban zones (Naachtun and Tikal). Fig. 8 Cluster analysis of density patterns. Ward linkage method, squared Euclidean distance. Such variety cannot be explained by sole reference to topographic, pedological, or ecological conditions; political and economic forces, including conflict, likely exerted commensurate or stronger influence. Furthermore, nearly all the areas identified as urban were also those incapable of sustaining themselves agriculturally, whereas all the rural and periurban areas were capable of producing food surpluses to feed the rest of the sample’s population. These data suggest that lowland Maya settlement was characterized by both high variability and high interdependence along a rural-urban continuum rather than by a dispersed, homogeneous, low-density settlement (94). A word of caution is warranted about assuming contemporaneity of all structures in this classification. For instance, three survey blocks—Zotz, Uaxactun, and Holmul—hold substantial centers that were abandoned in the Late Preclassic period and only partly resettled in the Classic period. However, because these blocks already form a statistically distinct group, we are confident that they would remain clustered even after factoring out all earlier settlement in those three areas. A few other patterns bear mentioning. First, the city of El Perú-Waka’ with a sustained density of 1100 structures/km2 (95)—nearly double that of Tikal—exhibits a uniquely nucleated urban core that transitions rapidly to a rural hinterland. The singular urban morphology of El Perú-Waka’, although it generally fits the urban area pattern, might reflect distinct economic or political influences on its development and organization. Second, although the Corona block was essentially rural, it did include several small pockets of urban density. This region had a strategic role in the expansion of the hegemonic Kaanul polity during the sixth through eighth centuries CE (96–102), so settlement nucleation here was probably conditioned by regional geopolitics in addition to local demographic processes. Finally, the urban and periurban settlement zones of Tikal extend over at least 76 km2, representing one of the largest continuous settlement zones in the Maya Lowlands (26).

Engineering and infrastructure Lidar data also elucidate the extent of lowland Maya investment in water management, regional communication, and defense. The Maya constructed reservoirs that required considerable labor (103). At the household level, reservoirs, quarries, wells, and underground cisterns were commonly cut into bedrock to collect rainwater (104–110). Within urban centers, the ancient Maya deepened, dammed, or bermed natural depressions to capture rainwater running off stucco-paved surfaces (111). Lidar data now show a greater extent of such large features. Large, round, and bermed reservoirs were built within wetlands to serve the needs of the rural population, and drainage ditches were cut into the edges of reservoirs to control flow during the wettest periods. Lidar-derived terrain models allow for a volumetric assessment of reservoirs, in measures likely to be conservative because of later infilling. Here, we present data only on the built reservoirs of large scale (i.e., those serving more than a single household) (29). In our sample, Tikal’s urban core had the largest amount of water storage capacity in human-made reservoirs. The Tikal Palace Reservoir alone could store 31,000 m3 of water, which would be sufficient to supply the inhabitants of the urban core for an entire year. Earlier cities such as Cival and Tintal were near natural depressions that provided generous amounts of water with little labor investment. For example, Tintal ringed a 2-km-wide sinkhole that could have contained more than 3 million m3 of water in wet years. A 2.5-km canal carried overflow into a natural drainage, preventing rising waters from flooding the city. In all cases, such features are monumental in scale and imply some form of centralized involvement in planning and execution. Notably, large-scale water infrastructure is not limited to city centers but occurs in the periurban and even rural zones surrounding major cities, and in some cases integrates cities and their peripheries. These findings indicate that the construction of monumental reservoirs and, less commonly, canals in both city centers and rural areas involved some form of centralized, institutional coordination. Yet this observation does not negate the possibility that smaller-scale infrastructure was independently built by households and corporate groups (104). Rather, the PLI data suggest that water management was neither fully centralized (112) nor left completely to individual households. Lidar data also reveal investment in regional and local interconnectivity. “Causeways,” elevated and paved roads of varied size, demonstrate economic and political integration by revealing formal connections between cities, smaller communities, and dispersed populations (113). Until recently, knowledge of formal political connections was restricted to epigraphic records of dynastic interaction and a small number of causeways (25, 26, 114–117). The lidar sample identified ~106 linear kilometers of causeway construction dating to the Late Preclassic and Classic periods. Intersite causeways are primarily associated with cities that rose to prominence during the Preclassic period (Tintal, Cival, San Bartolo) or Early Classic (Naachtun) and are surprisingly absent from the Late Classic period. They can reach up to 22 km in length (Tintal-Mirador) and 10 to 20 m in width, linking cities to smaller centers nearby. Examples of these occur at Tintal, where multiple monumental causeways radiate toward other centers. Intra-site causeways are associated with Late Classic period centers of all sizes. These roads are typically only a few hundred meters long, providing grand entryways to public spaces within a community. Tikal has the widest causeways (80 m) and the largest amount of internal paved causeway area (0.19 km2). Overall, causeways are more widely distributed than had been previously appreciated, especially in the Late Classic period. Across the sample, there is a general trend of increased causeway density from west to east and from south to north (Fig. 9), coinciding with regional gradients in settlement density, agricultural intensification, and water management infrastructure. Fig. 9 Map of survey blocks showing causeway density patterns. Causeways are shown in white; values represent cumulative causeway length (in meters) normalized by survey block area (km2). Perhaps the most obvious example of infrastructural investment comes in the form of defensive or military features pointing to a high incidence of conflict in the Maya Lowlands. Such fortifications are found at a scale and quantity (Fig. 10) matched only by the Tikal earthworks and the monumental fortifications at Becan, both known since the 1960s (118–120). To a notable extent, settlement density does not correlate with regional defense, in that some of the most densely settled blocks, such as Naachtun, have no defensive features. Individual defensive features—bridges, ditches, ramparts, stone walls, and terraces—were constructed as components of “built defensive systems.” These combined with natural defenses to protect “defended areas.” There were five types of built defensive systems: landscape ditch-and-rampart (type 1), hilltop ditch-and-rampart (type 2), contoured terrace (type 3), stand-alone rampart (type 4), and stone wall (type 5). These were unevenly distributed across the sample, reflecting variation in local geography, resources, and history. Fig. 10 Density of defensive features. Density is expressed as linear meters of defensive features per square kilometer normalized by survey block area. Built defensive systems are found in 36 discrete locations across the sample, with 31 meeting the criteria for defended areas (Table 10). In most cases (61.3%), defenses consisted of more than the establishment of perimeter walls. One-third (n = 10) of all defended areas displayed more defensive features within the perimeter than on the perimeter itself (Table 11). The Holmul and Zotz survey blocks exhibited the greatest effort in reinforcing defenses, suggesting a heavy investment in the protection of particular sites in those areas. We also identified “refuges,” highly defensible areas lacking any other perceptible architecture. Presumably, these were only used as a last resort and for brief stays. Clear examples appear at Kanalna and the El Achiotal promontory. Table 10 Defensive data for individual defended areas. Summary of defensive data for individual defended areas and perimeter statistics. Bold lines separate different lidar survey blocks (in order as Yala, Tikal, Uaxactun, Holmul, Zotz, Xultun, Perú, Corona, and Tintal). View this table: Table 11 Statistics for reinforcing defenses. Bold lines define natural breaks (Jenks) in the percentages of additional defenses in relation to an area’s perimeter. View this table: Several metrics were calculated to interpret the defensive data (29). We assessed the “perceived need” for defense of an area according to the percentage of the perimeter defended by built or natural defenses (Table 10), as well as how much defensibility factored into the siting of a settlement. The latter statistic is a preparedness index that factors out the natural defensibility of a location (Table 12). At some larger sites, such as Tikal and Tintal, substantial investment in building perimeter defenses around the urban core suggests that the need for such constructions emerged well after the initial establishment of the settlement. Finally, we assessed the cost of building defensive systems in relation to the total protected area (Table 13). The defensive features of greatest cost are located in the southern part of the sample, particularly the Early Classic Buenavista Valley linking El Zotz and Tikal. Table 12 Degree of preparedness for each defended area. The degree to which the need for defense was anticipated at the time of initial settlement. Scores closer to 0 indicate minimal preparedness. View this table: Table 13 Cost index for each defended area. Cost statistics for defending areas in relation to idealized forms. View this table: At the regional level, the Zotz, Tikal, Holmul, and Perú blocks show the greatest concentration of defensive features (Table 14). However, both the Zotz and Tikal regions demonstrate the largest investment of built and combined defenses in the sample. By controlling for survey block size, a line density analysis (Fig. 10) highlights the higher investment in defense along a west-to-east route. The trend continues to the east with several fortified hilltops in the Holmul region. The greatest density of high-cost features extends along a corridor associated with a Central Mexican (Teotihuacan) military intrusion into the Maya Lowlands in the late fourth century CE (121). Table 14 Regional defensive statistics by lidar survey block. TLD, total length defended; TBD, total built defenses. View this table:

Concluding remarks Since its first application in Maya archaeology (122), lidar technology has enabled the evaluation of lowland Maya society on a large scale. The PLI results confirm that lidar technology represents a watershed event in archaeological survey of forested environments. In the central Maya Lowlands, lidar survey will become indispensable for settlement research because of (i) the speed of the data-gathering process, (ii) the degree of detail attainable over large areas, and (iii) the ability to discern large and small features routinely undetected by traditional methods. The PLI lidar survey provides a uniquely large and continuous dataset for the central Maya Lowlands that is replete with evidence for ancient structures, canals, terraces, causeways, and defensive features. These data also show great variability in the intensity and distribution of these features over space. In sum, the PLI data unambiguously support the notion that the lowland Maya constructed a variable and contentious landscape in which a regionally interconnected network of densely populated and defended cities was sustained by an array of agricultural practices that optimized land productivity, resource diversity, and sustainability on a much grander scale than previously thought. These preliminary conclusions should encourage research that, among many possibilities, (i) quantifies labor investments and land productivity at every scale; (ii) addresses how projects of agricultural intensification were organized and managed, with implications for long-term sustainability; (iii) assesses the extent of economic interdependence within large regions by studying differences between urban and rural populations; (iv) determines the relationship among settlement densities, land use, and monumental infrastructure; and (v) develops models for the types of warfare implied by defensive features. Finally, by detecting most looters’ pits and assembling full inventories of ruins, lidar data offer a profound tool for cultural heritage management. Quantification of biomass and monitoring of biodiversity—an underexplored product of these lidar data—will also facilitate planning and management of regional conservation efforts and deepen our understanding of the relationship between ancient and modern environments.

Methods The Pacunam Lidar Initiative collected lidar data over an area of 2144 km2 of the Maya Biosphere Reserve (MBR) in northern Guatemala. These data were provided to the research consortium as bare-earth terrain models with grid resolution of 1 m. Given the coverage scale and resolution, these raster models make it possible to interpret aspects of ancient Maya society as warranted by visible remains of archaeological settlement. These data were analyzed to (i) estimate regional population levels, (ii) calculate agricultural capacity, (iii) determine settlement density patterns, and (iv) evaluate the extent of infrastructural investment in water storage, regional communication, and defense. Terrain data and its derivatives were mostly manipulated and analyzed through GIS applications (ArcMap 10.4 and 10.6, QGIS 2.18, and GRASS 7.2, and the Relief Visualization Toolkit 1.3), enhanced by supplemental calculations and statistical analyses performed in Microsoft Excel and SPSS. Terrain models were visualized using multiple methods. Archaeological features were digitized manually, using existing site maps as training samples, augmented by targeted ground-truthing during the 2017 field season. All PLI member projects uploaded their manually derived feature datasets to a shared cloud server, where they were collated to create a single regional dataset for each feature class (e.g., buildings, causeways, terraces, etc.). Population estimates are based on the total number of structures in the sample and include adjustments for Late Classic occupation, nonresidential structures, structures not visible on the surface, contemporaneity of occupation within the Late Classic period, and number of persons per structure. Settlement density classes were calculated using the Heatmap tool in QGIS 2.18 and classified according to published descriptions of settlement density in the cores and peripheries of numerous sites in the central Lowlands. Agricultural capacity was estimated by combining zones of intensive cultivation (defined on the basis of built agricultural features) with sparsely settled and unmodified uplands to define the total area available for agriculture; this area was then multiplied by a standardized estimate of traditional Maya swidden productivity. Water storage capacity was calculated using the r.lake function in GRASS 7.2. Causeways were classified using published typologies and statistics calculated by “calculate geometries” in ArcMap 10.4. Defensive features were classified by an original typology, their distribution was analyzed using the Line Density tool in ArcMap 10.6, and inferences regarding their cost effectiveness and other metrics were derived from statistical analyses performed in Microsoft Excel. See (29) for detailed methodology and analytical procedures for each of these analyses.

Supplementary Materials www.sciencemag.org/content/361/6409/eaau0137/suppl/DC1 Materials and Methods Supplementary Text Figs. S1 and S2 Tables S1 to S4 References (123–144)

http://www.sciencemag.org/about/science-licenses-journal-article-reuse This is an article distributed under the terms of the Science Journals Default License.