A Low-Cost GPS-Based Protocol to Create High-Resolution Digital Elevation Models for Remote Mountain Areas

Abstract Researchers and development practitioners in remote mountain areas rely on elevation data to study vegetation dynamics, natural hazards, land use, and other patterns. However, despite advances in technology, accurate digital elevation models (DEMs) with spatial resolution <30 m do not exist for most of the world's montane regions. We used a low-cost GPS-based protocol to construct a high-resolution (10 m) DEM for a rugged, remote mountain site in the northern Peruvian Andes. Elevation data were collected with handheld GPS units and combined with digitized and interpolated points within a Geographic Information System to generate a 10 m DEM. Additional DEMs were generated using 50%, 20%, and 15% of the surface points collected and from a 1∶100,000 topographic map and ASTER GDEMv2 data. Estimated absolute vertical accuracy of the GPS surface-point DEMs was significantly lower than that of the ASTER GDEMv2 and topographic map DEMs. Relative vertical accuracy, a better measure of DEM quality, was considerably lower for all 6 DEMs than absolute vertical accuracy. Depending on project budget, time, and labor availability, this method can be used to produce DEMs with high spatial resolution and substantially improved relief maps for research, visualization, and communication purposes. Implementation of this method is practical in locations without access to electricity or post-processing correction facilities, open-canopy land covers, and projects with small budgets that involve local participants.


Introduction
Coupled with latitude and longitude, elevation provides 3-dimensional (3D) locational information describing terrain which is essential in mountain research and development (Kö rner 2007;Malhi et al 2010). Alexander von Humboldt was one of the first Western explorers to recognize this: his expeditions in montane Mexico, Colombia, and Ecuador demonstrated that knowing the 3D location of biophysical features on the Earth's surface was important for cartography and for understanding distributional relationships among biotic, abiotic, and human factors interacting along elevational gradients (Godlewska 1999;Zimmerer 2006;von Humboldt 2013). Anthropologists, geographers, and ecologists have since sought to quantify and visualize how elevation influences diverse phenomena in montane regions (McVicar and Kö rner 2013). As examples, studies show that increasing elevation causes fundamental changes in species distribution (Feeley et al 2011), crop diversity (Zimmerer 1999), agricultural land use (Guillet 1981;Brush 1982;Young 1993a), net primary productivity (Beck et al 2008;Zhang et al 2013), and biogeochemical cycling (Girardin et al 2010;Ramos-Scharró n et al 2012). Moreover, elevation and topographic position affect the vulnerability of mountain communities to natural hazards, disease, food insecurity, and climate change (Larsen and Torres-Sá nchez 1998;Huddleston et al 2003;Young and Lipton 2006;Siraj et al 2014).
Current understanding of elevation, and more broadly of topography, as a determinant of landscape pattern and process in montane regions has improved considerably over the past 3 decades due in large part to the advent of digital elevation models (DEMs). Freely available DEM data now enable numerous terrain parameters (eg slope and aspect) to be derived for geospatial analysis and modeling (Pike 2000). For example, the Shuttle Radar Topography Mission (SRTM, www2.jpl.nasa.gov/srtm/) provides 90 m resolution elevation data ( Figure 1A) suitable for regional and continental applications, such as mapping land cover and interpolating climate surfaces (Hijmans et al 2005). Also freely available, the Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model Version 2 (ASTER GDEMv2, http://asterweb.jpl.nasa.gov/gdem.asp) is a data set with 30 m intervals ( Figure 1B) that has been  (Saadat et al 2008;Khan et al 2014) and to investigate geomorphic and hydrologic processes at landscape to regional scales (de Vente et al 2009;Santini et al 2009). Both SRTM and ASTER GDEMv2 are widely used and provide near-global coverage. However, SRTM contains no-data or void areas, which must be filled before processing, and ASTER GDEMv2 contains anomalous values and artifacts (see Figure 1B) that must be filtered or modified using interpolation methods (Santini et al 2009).
Newer technologies such as light detection and ranging (LiDAR) allow creation of submeter resolution DEMs ( Figure 1D). Nonetheless, free or low-cost data from LiDAR sensors are restricted to a few locations, and data acquisition may be cost prohibitive for most researchers and development practitioners working in remote mountain areas. Unmanned aerial vehicle (UAV) photogrammetry is practical for areas covering less than a few hectares, such as archaeological sites (Lambers et al 2007;Verhoeven 2009), although in many areas their use remains impractical because of platform, sensor, operating, and environmental constraints (Anderson and Gaston 2013). Photogrammetry from high-resolution sensors, like Digital Globe's World View, GeoEye satellites, and the Advanced Land Observation Satellite (ALOS) PRISM (Panchromatic Remote-Sensing Instrument for Stereo Mapping), is another possible source of ,5 m DEMs. Unfortunately, atmospheric conditions in tropical and high mountain regions coupled with steep terrain often preclude use of this technique to develop suitable DEM data (Kä ä b 2002). In sum, although capabilities for creating DEMs are changing rapidly, grids with spatial resolution ,30 m simply do not exist for many remote montane regions around the world.
This gap in DEM spatial resolution ( Figure 1C) is problematic because several issues and priorities in mountain research and development require better understanding of processes influenced by topography at spatial scales ,30 m. For example, higher-resolution data are critical to predict landslide hazard (Stark and Hovius 2001), a widespread phenomenon in humid highlands partly controlled by slope angle. A study in Puerto Rico found that average landslide width was 14 m, less than half the size of one ASTER GDEMv2 pixel (Larsen and Santiago Romá n 2001).
Alpine plant response to climate change at high elevations may also vary over small spatial scales (#10 m) as a result of terrain roughness and gradient (Young 1993b;Coblentz and Keating 2008). As well, fine-scale topographic data are needed to quantify water storage and design water control structures in environmental services programs (Guzha and Shukla 2012) and to assess sediment, water, and nutrient fluxes in hillside agroecosystems (Ericksen et al 2002). Unfortunately, the coarse spatial scale of freely available DEMs cannot be used in such assessments.
In addition to their use in research, DEMs can be employed to render 3D representations of surface topography that are useful in a wide variety of applications, including urban and hydrologic modeling; glacier, landform, and soil classification; and river bank and forest management (Liu 2008). Limited availability of high-resolution data affects the scale at which these elevation-dependent phenomena can be examined, visualized, and applied.
Given these challenges, the objectives of this study were to (1) develop a high-resolution (#10 m) DEM for a remote site in the northern Peruvian Andes using a lowcost protocol based on the Global Positioning System (GPS); (2) compare and contrast the spatial resolution and accuracy of DEMs generated from surface-point, contourline, and remotely sensed data; and (3) identify advantages and disadvantages of a GPS-based protocol for constructing high-resolution DEMs in remote mountain areas, thus making these findings of general importance.

Study area
This study was conducted in Río Abiseo National Park along the eastern slopes of the Andean Eastern Cordillera in northern Peru ( Figure 1A-C). The park encompasses 274,520 ha along a vegetational and altitudinal gradient spanning 700-4200 m. Three ecological zones exist within the park: tropical premontane forest (800-1500 m); tropical montane forest (1500-3500 m); and tropical alpine grassland above timberline (3500-4500 m) (Young 1993b). The tropical alpine zone is separated from closed tropical montane forest by a relatively narrow timberline belt in which numerous small (,5 ha) forest patches are distributed (Young 1993a(Young , 1993b; Figure 2A, B).
This protocol was implemented and evaluated in Callejó n Rojas (7u56.8879S, 77u21.3579W), a north-south trending U-shaped valley (670 ha) located in the timberline zone ( Figure 2) and ranging from 3390 to 4070 m. Long-term climate data are not available for this site. However, remotely sensed precipitation data (CPC MORPhing technique) for 2003-2008 indicate an approximate annual rainfall of 1230-1630 mm, with 80%-90% of rain falling between October and April and a drier period between May and September. This site is also exposed to orographic fog throughout the year (Young 1993a(Young , 1993b. As is typical of remote mountain sites, high-resolution topographic mapping in Callejó n Rojas involves several logistical and methodologic challenges. First, few cloudfree satellite images are available for this region. Second, recent topographic maps do not exist; 1:100,000 maps produced by the Instituto Geográ fico Militar de Perú from 1962 aerial photographs are the only topographic maps available ( Figure 1C). Third, total travel time on foot from the valley to the nearest source of electricity (a park guard station) is approximately 6 hours. The valley itself is inaccessible to vehicles, thus using bulky survey and field-mapping equipment is not feasible. Fourth, lack of electricity limits the use of rechargeable batterypowered equipment such as a more advanced GPS

Methods
Elevation data used to produce DEMs can be remotely sensed or derived from contours or surface points (Hutchinson and Gallant 1999). All 3 types of data were used in this study.

Elevation and landscape data collection
From 10 to 19 July 2011, handheld GPS units and a camera were used to collect elevation and ancillary data for the construction of a high-resolution DEM of Callejó n Rojas. Local residents with knowledge of the area were hired to assist with GPS data collection (Figure 2A, C, D). Images of the valley, including an ALOS image with ,2.9 m resolution (16 May 2007) and aerial photographs with ,2.5 m resolution (26 December 1962) were used during initial ground surveys to familiarize the team with the geography of the valley (Figure 2A). Field assistants were trained in the use of GPS technology, including basic equipment operation, best practices for recording data, and an overview of satellite paths and interference. Because of the lack of electricity and postprocessing capabilities, standard consumer-grade, AA-batterypowered, handheld GPS units-a Garmin GPSMAP 62st and 2 Garmin eTrex H units-were used to collect elevation points.
Because the density and distribution of point data can affect DEM accuracy (Chaplot et al 2006;Weng 2006;Erdogan 2009), GPS points were spatially distributed across the valley. First, 13 subwatersheds were delineated using the ALOS image, aerial photographs, and field observations. Second, specific areas were assigned to each of the surveyors for GPS point collection. Third, GPS units were set to 3-second track intervals, allowing multiple points to be taken in rugged terrain. GPS units recorded track points with horizontal and vertical location information ( Figure 3A). Fourth, where possible, track points were collected parallel and perpendicular to slope to obtain a broad distribution of elevation points. Data collection ceased during rain and fog events.
In addition to GPS track points (n 5 33,138), ground control points (GCPs) (n 5 365) were registered at locations (eg peaks) and features (eg lake edges) discernible in aerial photographs and representative of distinct land-cover types ( Figure 3B). Multiple photographs were taken at each GCP location with a Canon Power Shot SX120 IS. Land cover, photograph direction, and camera angle were recorded. To enhance positional accuracy, all GCPs had an uninterrupted GPS signal from at least 4 satellites and an average horizontal error reading of ,3.5 m.

DEM construction and accuracy validation
GPS track points, GCPs, and photographs were downloaded onto a desktop computer. A Geographic Information System (GIS) (ESRI 2008(ESRI , 2012) was used to visualize, digitize, edit, interpolate, and analyze spatial data. Track points and GCPs were merged into a single point shapefile ( Figure 3A). A 5 m buffer was created around each GPS point and a random point shapefile was generated with 34,000 points to fill spaces between buffers. Using the Editor tool in ArcMap, elevation values were assigned to each random point using a manual interpolation method. Track point data, field notes, photographs taken from multiple vantage points, aerial photographs, and the 2007 ALOS image were used to estimate the elevation of each random point. Specifically, elevation at each point was interpolated using adjacent elevation points, with closer points weighted more heavily. Along slopes, we used average elevation between upslope and downslope points to assign point elevation values. Field photographs from higher elevations provided oblique aerial perspectives of the valley that proved particularly useful for gauging changes in relief.
To resolve errors in vertical position due to satellite geometry and ionospheric effects, GPS points were reviewed and edited during different stages of the GIS work. For example, our first step in the GIS stage was to assess the track point data. Points were interpolated using ordinary kriging, which averages subsets of neighboring points to derive values for unmeasured locations (ESRI 2012), to identify points with significant vertical error. These points appeared as deep pits or peaks in the interpolated model and were removed from the data set. During the digitizing phase, points with elevation values markedly different from those of surrounding points were edited to reflect topography as it appeared in field photographs and/or was recorded in surrounding points. When different GPS units recorded similar elevation for points in close proximity, it was assumed that vertical and horizontal position error was low. After all elevation points were edited and assigned using the manual interpolation method, ordinary kriging was again used to interpolate point data and to further identify and edit false pits and peaks.
GPS track points, GCPs, and manually interpolated points with elevation data were merged into a shapefile containing 66,877 elevation points. GCPs (n 5 365) were not included as these were used for DEM accuracy validation. Elevation data were used to produce a DEM with an output cell size of 10.0241 m. Ordinary kriging was used because interpolation methods have been found to result in little difference among DEMs at micro to catchment scales when sampling densities are high (Chaplot et al 2006). Three additional DEMs were generated to assess how differences in time allocated to field and computational work would have affected DEM resolution and accuracy. For these DEMs, 50%, 20%, and 15% of the data points were randomly selected and used in DEM construction. The same interpolation approach was used to allow direct comparison of DEMs. The spatial resolution (m) of each DEM was determined using Tobler's (2000) spatial theorem equation: where d 5 the dimension of the region of interest.
DEMs of the valley were produced using ASTER GDEMv2 and a 1:100,000 digitized and rectified topographic map with 50 m contours (Instituto Geográ fico Nacional de Perú 1984). ASTER GDEMv2 has a cell size of 30 m ( Figure 1B). For the 1:100,000 topographic map ( Figure 1C), contour lines were digitized and a triangulated irregular network was created, from which a DEM was then generated. The output cell size for this contour DEM was 38.1 m.
We compared the accuracy of the 6 resulting DEMs using 3 methods. The difference between recorded and modeled elevation at each GCP was calculated, and the distribution of differences was examined with histograms. The absolute vertical accuracy of the 6 DEMs was computed using a standard measure of map accuracy, the root mean square error (RMSE) (Fisher and Tate 2006): where z DEM is DEM elevation, z REF is GCP elevation, and n is the number of GCPs. We also calculated the relative, or point-to-point, vertical accuracy of the DEMs (Gesch 2007). The point nearest to each of the 365 GCPs was identified. The relative vertical accuracy (RV) for each point pair was then calculated using the following equation (National Digital Elevation Program 2004): where D REF is the elevation difference between GCP point pairs and D DEM is the elevation difference between the modeled point pairs.

Results and discussion
We developed 4 DEMs with ,30 m spatial resolution using a low-cost GPS-based protocol for a remote mountain site in the northern Peruvian Andes, for which elevation data at this resolution do not currently exist. Three-dimensional visualizations of the GPS surfacepoint DEMs (hereafter GPS DEMs) show substantial improvements over the ASTER GDEMv2 and the contour DEM for this study site (Figure 4). The ASTER GDEMv2 and contour DEM underestimate elevation more than the GPS DEMs and have vertical deviations that are 5 to 9 times larger than the GPS DEMs ( Figure 5). For example, the contour DEM overestimates elevation by up to 185 m, whereas the GPS DEMs overestimate elevation by no more than 31 m. With regard to absolute vertical accuracy, the  (Table 1). For this site in northern Peru, the RMSE of the ASTER GDEMv2 is greater than the reported average global accuracy of 620 m and similar to the RMSE estimated for a site in the Andes of Argentina (630 m) (Eckert et al 2005). The high RMSE of the ASTER GDEMv2 we report is likely due to the steep topography of our study site, inaccuracy of our GCP points, or a combination of these factors. Our use of simple handheld GPS units prevents us from knowing the absolute vertical accuracy of the GCPs. Also, several studies have shown that the RMSE of the ASTER GDEMv2 data is positively correlated with slope angle (Toutin 2002;Fujita et al 2008;Wang et al 2012).
The GPS DEMs generated for Callejó n Rojas have much lower RMSEs than the ASTER GDEMv2 and contour DEM because GCPs were used as part of the manual interpolation process. However, we cannot determine the absolute accuracy of the GPS DEMs given the unknown vertical error associated with track points and GCPs. Lack of access to electricity, cost associated with signal correction services, and weight constraints limited our ability to use differential GPS. In addition, this region falls outside the North American Wide Area Augmentation System (WAAS) broadcast area. WAAS provides GPS signal corrections that can improve positional accuracy by ,3 m. Without WAAS correction, consumer-grade GPS units have a 95% confidence interval of ,15 m for horizontal and ,22.5 m for vertical position (Pawlowiciz 2007). Position measurements can also vary over time because GPS satellite geometry is constantly changing (Pawlowiciz 2007).
Interpolation algorithms used to generate DEMs represent another source of error in the modeled elevation values (Fisher and Tate 2006). Ordinary kriging is a standard interpolation method (Yilmaz 2007) but may not be the most appropriate when sampling densities are on the low side (Chaplot et al 2006). Despite these caveats, the magnitude of difference in estimated absolute vertical accuracy among the DEMs underscores the limitation of the ASTER GDEMv2 and contour DEM for small-scale research in steep terrain. Nonetheless, both approaches are widely used in remote mountain regions as they represent some of the only sources of elevation data available (Chirico et al 2012).
For researchers and development practitioners working in remote montane regions, relative vertical accuracy-the accuracy with which a DEM represents relief-may be a more useful metric of DEM quality than absolute vertical accuracy. This is especially true for terrain features such as slope and aspect, which are derived from DEMs and calculated as the difference between adjacent elevation values (Gesch 2007). All the DEMs of Callejó n Rojas have relative vertical accuracies that are improved over absolute accuracies (Table 1). This indicates that GPS surface-point DEMs can be used to create fine-scale maps of topographic features, such as slope and aspect, and to calculate indices, such as topographic wetness index, that rely on improved estimates of local relief. Such maps can also be used for visualization and communication purposes. The coarserscale DEMs similarly show considerable improvements in relative vertical accuracy, although at a much lower spatial resolution (Table 1).
The method we propose here has several obvious limitations. Potential sources of error include insufficient accuracy, density, and distribution of source data (Gong et al 2000). As discussed, handheld GPS units have horizontal and vertical errors. Variable topography and high relief in mountain environments affect satellite signals as well. In addition, we did not collect track points or GCPs along cliffs, extremely rocky terrain, or other unsafe areas. For this reason, minimizing human error during point collection and manual interpolation, as well as collection of detailed notes, are critical for developing accurate DEMs using this method. Depending on the desired resolution, implementation of this protocol can be labor intensive, both in the field and in the lab. The 10 m resolution GPS DEM involved collection and manual interpolation of 66,877 points and required a total of ,644 person-hours, not including training and travel time (Table 1). Moreover, GIS software is a prerequisite for this protocol as it is for generation of DEMs from remotely sensed and contour-line data. For some GIS applications, open-source mapping software may be an alternative to ArcMap, whereas Google Earth and Bing Maps offer free satellite imagery and some capabilities for producing 3D visualizations.
Our analysis also highlights several advantages of the GPS-based protocol for remote mountain regions. Assuming a computer and mapping software are available, high-resolution DEMs can be constructed at a relatively low cost compared with other methods, such as LiDAR and UAV photogrammetry. The GPS units, batteries, supplies, and camera cost approximately US$890. With respect to labor, the highest-resolution (10 m) GPS DEM required walking more than 200 km for 144 person-hours (Table 1). Assuming an 8-hour workday, 4 people could collect these data in ,4.5 days.
The primary advantage of our protocol is that researchers and development practitioners can weigh spatial resolution requirements against project budget, availability of personnel, and time. For example, a 50% reduction in effort (ie collection and interpolation of 33,439 surface points) results in only a 4 m decrease in spatial resolution and ,1 m decrease in relative vertical accuracy compared to the 10 m GPS DEM (Table 1). Yet, the ,14 m GPS DEM would still involve a total of 322 person-hours. An 80% reduction in effort results in the generation of a much lower-resolution (,22 m) DEM, but this surface-point DEM is still 8-16 m higher in spatial resolution than the ASTER GDEMv2 and the contour DEM. In addition, the ,22 m GPS DEM has a much better relative vertical accuracy than the remotely sensed and contour DEMs. Use of only 15% of the surface points causes a disproportionate decrease in spatial resolution. In sum, the 20% effort scenario appears reasonable in situations where there are existing budgetary or time constraints. Another advantage of the GPS method over the ASTER GDEMv2 and contour DEM is that it provides opportunities for local participation in the research process ( Figure 2). In local communities, outreach and technology transfer actively involve people in science and can foster community empowerment (Dana 1998;Offen 2003). With inexpensive GPS units now widely available and GIS increasingly available, proper training in the use of geospatial technology and manipulation of geospatial data is particularly important in remote mountain regions (Bussink 2003;Heinimann et al 2003).
Through long-term fieldwork, researchers and development practitioners become intimately acquainted with field sites yet often lack access to high-resolution DEMs. As a result, they must rely on low-spatial-resolution DEMs even at sites where these may not be especially useful. We demonstrate how intimate knowledge of a study site coupled with elevation data can be integrated within a GIS framework to produce a high-resolution DEM for places for which such data are currently not available.

Conclusions
Few high-resolution DEMs exist for tropical and other high mountain regions, resulting in limited data accessibility for mountain researchers and development practitioners. DEMs can be constructed using surfacepoint, contour-line, or remotely sensed data. In this study, GPS collection of elevation points augmented with GIS manual interpolation of elevation based on detailed field notes, photographs, and ancillary data were used to generate a 10 m resolution DEM of a remote Andean site in northern Peru. DEMs produced using this method have absolute elevation accuracies that are subject to uncertainty; however, relative vertical accuracies suggest that this is a method with significant utility for obtaining terrain parameters such as slope and aspect that rely on differences between elevation points. DEMs produced using this method also cost substantially less than DEMs developed from LiDAR or UAVs, and by involving local participants in data collection, GPS-based DEMs can foster outreach and technology transfer. Our final product, though not free of error, provides a substantially improved relief map and represents a type of solution that may be generally useful in many other situations.

ACKNOWLEDGMENTS
We are grateful to Don Pancho, Elias, José, Javier, Carlos, and Eber for their assistance with field data collection. To all the park guards of Río Abiseo National Park we owe special thanks for ensuring our safety in the field and during travel. Thanks to Susy Castillo, Blanca León, and Asunción Cano for their tireless assistance with the coordination of this project. This research could not have been conducted without the support of SERNANP (permits Nu 001-2009-SERNANP-DGANP-PNRA, Nu 002-2011-SERNANP-PNRA), Río Abiseo National Park staff, the Museo de Historia Natural in Lima, and the Asociación Peruana para la Conservación de la Naturaleza (APECO). We are eternally indebted to Mariella Leo for assisting us with coordination and logistics. This research was funded by the National Science Foundation Minority Postdoctoral Research Fellowship (SBE #0905699, to A.G. Ponette-González), the National Geographic Society Committee for Research and Exploration (#8877-11, to A.G. Ponette-González), and the United States National Science Foundation (NSF) (DEB #1146446, to K.R. Young).