INTRODUCTION

High resolution topography is needed in areas with escarpments and fractures caused by different subsidence and uplift rates, as is the case for a 1.44 km² maar lake bottom located in Parangueo, Mexico, where traditional surveying methods can not provide the required detail to characterize all topographic features. A possible solution to this problem would be the generation of a LiDAR point cloud, but its high cost hinders its application. In recent years, the acquisition of high resolution topography data has become easier than before due to two new tools: Structure from Motion algorithms and light Unmanned Aerial Vehicles (UAVs). Structure from Motion (SfM), is a computer vision technique that involves the simultaneous recovery of 3D camera motion and 3D scene structure from a collection of tracked 2D features on overlapping pictures (^{Szeliski, 2010}). When applied to aerial photography, this technique provides data quality and resolutions comparable to-or better-than those obtained with LiDAR and classic photogrammetry with unprecedented ease of use and low cost (^{Fonstad et al., 2013}; ^{Leberl et al., 2010}). Furthermore, SfM estimates both internal and external camera orientation parameters, including nonlinear radial distortions, making it possible to use consumer grade digital cameras to develop 3D point clouds. These cameras can be carried by Unmanned Aerial Vehicles (UAVs), which are classified as fixed wings or multirotors and can be flown in manual or programmed modes. Because of its potential, the combination of these two technologies has been explored on a wide range of studies, including archeology (^{Verhoeven et al., 2009}), rangeland assessment (^{Laliberte, 2009}), vegetation structure (^{Dandois and Ellis, 2010}), biodiversity in forests (^{Getzin et al., 2011}), weed management (^{Torres-Sánchez et al., 2013}), ecohydrology (^{Vivoni et al., 2014}), precision agriculture (^{Verger et al., 2014}), glacier dynamics (^{Immerzeel et al., 2014}; ^{Ryan et al., 2015}) and landslide deformation and erosion (^{Stumpf et al., 2015}).

Because these are relatively new technologies, different works have been conducted to assess their accuracy and/or to identify which of the available SfM softwares yield better results, as listed in Table 1: ^{Johnson et al. (2014)} and ^{Ouédraogo et al. (2014)} compared the point cloud obtained with SfM to that of a terrestrial laser scanner, while ^{Ouédraogo et al. (2014)}, ^{Sona et al. (2014)} and ^{Turner et al. (2014)} compared the results obtained from different softwares that create 3D point clouds using SfM algorithms, such as PhotoScan, MicMac, Pix4D and Bundler. In order to undertake these studies, digital cameras have been mounted on hellium balloons, fixed wings, and multirotors, with flight elevations above ground that range from 30 to 500 meters (Table 1); the one study that was flown above 500 meters used a tripulated helicopter. Of the 16 studies listed, the one by ^{Ryan et al. (2015)} covered the largest area (5 km²⁾ using a fixed wing to study glaciar calving dynamics in Greenland. However, their study is also the one with the largest georeferencing errors, with RMSE values of 141 cm and 190 cm on the horizontal and vertical; the large errors of their study were caused by a lack of measured Ground Control Points (GCPs), because they used identifiable geographical features on orthophotos as GCPs. On the same table, it can also be seen that multirotors have been used to cover small areas, as the largest area covered so far with a multirotor was of 0.275 km² by ^{Tonkin et al. (2014)}, who developed a Digital Surface Model of morraines and who reported horizontal and vertical RMSE values of 3.3 cm and 1.8 cm; of note is that they used a mirrorles camera (Canon EOS-M), with an image sensor identical to that used on traditional DSLR Canon cameras (22.3 mm × 14.9 mm) such as the Canon 550D used by ^{Lucieer et al. (2014)}, Turner et al. (2014) and ^{Mancini et al. (2013)} but lighter (0.21 kg vs 0.53 kg without lens). However, ^{Sona et al. (2014)} used a point-and-shoot camera (sensor size = 6.20×4.65 mm) reporting acceptable RSME values both horizontally and vertically (1.0 cm and 2.0 cm, respectively).

Table 1 Topography-related studies where Unmanned Aerial Vehicles have been used. Two or more values on the error column reflect either errors for diferent flight dates, flight heights or Structure from Motion packeges used. 

^m The authors reported mean values, not RMSE values *The RMSE values are arranged according to the software packages listed in the previous column

From the studies developed so far (Table 1) it seemed that the use of Structure from Motion applied to images acquired by an UAV could be the solution to our task: the development of a detailed topographic map showing all structures casued by both uplift and subsidence in Parangueo. For our study area, the use of a fixed wing was discarded because there is not an adequate landing place; additionaly, we needed a multirotor that could be carried around the crater easily, as the area to be mapped was large compared to those areas previously studied with a multirotor (Table 1). The solution was the use of a small quadcopter together with a point-and shoot camera and the use of the PhotoScan SfM software to create a high resolution DSM. We also assessed the accuracy of the generated DSM, as the high albedo of the sediments could represent a problem on the image acquisition phase; the details of our study follow.

STUDY AREA

Rincón de Parangueo is a Quaternary maar located in the northcentral part of the Trans Mexican Volcanic Belt (^{Aranda-Gómez et al., 2013}), on the central part of Mexico (Figure 1). A total of four crater lakes (with Parangueo being the largest in extension) were found in the Valle de Santiago (Figure 1) and are now dry due to the large groundwater extraction in the irrigated croplands that surround the region. The gradual desiccation of the Parangueo lake between 1986-1995 can be seen on the false color composites of Figure 2, which highlights the irrigated areas in the region. The deformation inside Rincón is remarkable as it illustrates several important processes associated with gravity deformation, and the lake basin may be seen as a mesoscale model of gravity deformation in passive margins. The dry bed lake in Rincón de Parangueo displays a large number and variety of structures associated with active deformation (^{Aranda-Gómez et al., 2014}) caused by the large subsidence rate at the bottom of the maar, estimated at 50 cm/year (^{Aranda-Gómez et al., 2013}). Some of the structures inside Rincón de Parangueo are similar in form and origin to structures documented at the bottom of the Gulf of Mexico and in adjacent areas in Texas (^{Rowan et al., 1999}). Near the coast of the desiccated lake there is clear evidence of radial extension, whereas closer to the center of the lake the evidence is of radial compression (^{Cobbold and Szatmari, 1991}). A remarkable feature of the unconfined deformation in Rincón de Parangueo is that it is occurring on dry mud in contact with the atmosphere, deforming in a cataclastic fashion and gliding in nearly rigid blocks (^{Schultz-Ela, 2001}), with plastic deformation occurring below the surface, where the lake sediments are still wet. The most compelling evidence of both plastic deformation and spreading (^{Schultz-Ela, 2001}) in the system is the formation of mud-injection domes and folds in the area between the toe of the inner scarp and the center of the basin. Due to the relatively small size of the bottom of the crater (1.2 km in diameter), the development of a topographic map to show all the structures associated with its deformation has been a difficult task due to the high resolution needed to characterize all the aforementioned structures.

Figure 1 Location of the study area. The Parangueo Maar is located in the Valle de Santiago, in the central part of Mexico, where other maars were found before heavy groundwater extraction started in the mid 1980s; currently, only two of them have a small amount of water. The background image is a true-color SPOT-5 composite acquired on March 5th, 2007 and highlights the irrigated areas around the study area. 

Figure 2 Drying of the Parangueo lake as observed from three false color composites of Landsat-5 images acquired before the rainy season from 1986-1995: (a) April 1986, (b) March 1990, and (c) May 1995. The false color image (bands 7,4,2) highlights the irrigated areas near the study area. Coordinates are in meters, UTM-14. 

Topographic maps of Parangueo

When we started studying the deformation at the bottom of Rincón de Parangueo, the highest resolution data available was the LiDAR data released by INEGI, with a resolution of 5 metres; however, as shown on Figure 3, this resolution is not sufficient to map the escarpment, fractures and domes that appear at the bottom of the dry bed lake. To overcome this lack of information we attempted to develop a detailed topographic map from a survey of nearly 5000 points distributed at the bottom of the lake; however, these surveyed points where not enough to register the complex system of faults and fractures that appear in the lake sediments. Other efforts to document the topography included the use of ALOS Prism and Envisat images; however, these efforts were also fruitless, as the resolution of the images and the high albedo of the sediments hinders their application to this problem. The solution to map all the features of Parangueo was the use of a small quadcopter that carried a point and shoot camera as described in the following sections.

Figure 3 Lidar and SPOT data for Parangueo: (a) Shaded relief (sun azimuth=225°, elevation=30°) of a 5 meter resolution LiDAR obtained from INEGI and (b) True color composite of SPOT-5 imagery (2.5 m resolution, acquired on 2007-03-05, during the dry season) showing the high albedo of the dry lake bed sediments. Some terrain aspects start to appear with the LiDAR data, such as the limit of the lake and the platform; however, no analysis can be undertaken at the crater's bottom. Coordinates are in meters, UTM-14. 

METHODOLOGY

For this work, a Phantom 2 quadcopter (manufactured by DJI Innovations: http://www.dji.com/product/phantom-2) was used. Without propellers, the diagonal length of this quadcopter is 39 cm; with propellers, this length increases by 20 cm and can carry a point and shoot camera in addition to its battery, which has a capacity of 5.2 Ah, allowing the quadcopter to fly for about 13 minutes when carrying the setup used in this work. This setup is shown in Figure 4 and included: (a) a GPS watch (Garmin Forerunner 305) attached to the upper part of the quadcopter, (b) a Canon S110 camera pointing downwards and attached on the lower part of the quadcopter's body, and (c) DJI's 2.4 GHz datalink (installed inside the quadcopter's body) that allows communication between a PC and the quadcopter for flight programming. The Phantom 2 quadcopter has a GPS enabled flight controller (Naza-M V2) with accuracies of 2.5 m on the horizontal and 0.8 m on the vertical. The GPS watch is used to extract flight details (i.e. coordinates, elevation and speed) after all flights are completed; these extracted coordinates are then written on the picture's EXIF header by matching both the pictures and the coordinates acquisition time through the use of JetPhoto Studio. A Canon S110 camera was selected because it can be set in intervalometer mode by using the Canon Hacker's Development Kit (CHDK, http://chdk.wikia.com/wiki/CHDK). The CHDK is written on the SD card of the camera and is automatically loaded when the camera is turned on, allowing extra features not available on the camera's default firmware. Through the use of the CHDK, the camera was set to take pictures every two seconds in raw format, focused at infinity, ISO 500, a shutter speed of 1/600 seconds, a focal length of 11 mm (equivalent to 50 mm in a full frame camera, due to the S110's crop factor of 4.62) and an aperture of f/8 (the maximum possible on this camera); the SD cards used on the camera had a 32 GB capacity. The flights were programmed at an elevation of 1900 m a.s.l. because the bottom of the crater has an elevation that varies from 1686 to 1712 m a.s.l., while the elevation of the launching spots at the crater's rim varies from 1790 to 1850 m a.s.l. In order to guarantee a minimum side overlap of 60%, the flight lines had to consider the ground coverage of each picture, which requires the camera's sensor and pixel sizes (7.6 mm × 5.7 mm and 0.0019 mm, respectively). With this information, and assuming a flight elevation above ground of 200 meters, the Ground Sample Distance (GSD) can be by:

where PS is the sensor's pixel size, Fh flight height and Fl focal length. According to eq. 1, at a height of 200 m above ground, the footprint of each picture is a rectangle with sides of 136 and 102 meters. The software used in this work (Photoscan) requires a minimum overlap of 80% on the forward direction and 60% on the lateral, thus, the maximum distance between each flight line has to be 52 meters. Accordingly, the flight lines were set-up at 1,900 m a.s.l., a distance of 50 meters between them and a horizontal speed of 12 m/s (43.2 km/h); a total of 22 flight lines were required to cover the 1.44 km² crater's bottom as shown in Figure 5. These 22 flight lines were covered through 4 different launches, which required a battery replacement; in order to maximize battery life, the launching and landing points were not the same, as clearly illustrated on the second launch (Figure 5 and Figure 6). On the third launch, the battery started to run out of power and the quadcopter descended smoothly (Figure 5) as the flight controller is programmed to land when the battery reaches 25% of its capacity. Fortunately, we were able to find the quadcopter and continue with the flights as no damage occurred neither to the quadcopter nor the camera (even though it landed on a branch); accordingly, the fourth launch was set up so its last flight line would cover the final part of flight line 18 (last flight line of the third launch), as shown on Figure 5.

Figure 4 Set-up used in this study with a Phantom-2 quadcopter (diagonal length of 59 cm including propellers): (a) a portable GPS (Garmin Forerunner 305) attached on the top of the quadcopter in order to record flight data, (b) a Canon S110 placed on the bottom of the quadcopter pointing downward, and (c) DJI's ground station software was run on a portable computer in order to program all required flights-picture taken by the quadcopter during take-off. 

Figure 5 A total of 22 flight lines were required to cover the crater's 1.2 km × 1.2 km bottom. The flight lines shown were registered by the GPS watch attached to the quadcopter and visualized on Google Earth: (a) plan view, where the distance between flight lines can clearly be seen, and (b) perspective view showing both take-o↵ and landing points as well as the elevation of all flight lines; flight elevations were set at 1900 m a.s.l. The turns at the end of each flight line were set as adaptive bank turns and a total of four flights were required; the quadcopter was landed at the end of each flight to replace the battery. 

Figure 6 Details of launch 2 (flight lines 7-12, shown in green color on Figure 5) as recorded by the GPS-watch: (a) velocity, and (b) elevation. The quadcopter was programmed to use adaptive turn banking; with this configuration, the quadcopter does not stop on every turn, but reduces its velocity as clearly seen on (a), while maintaining its elevation (b), where the nearly 50 m in elevation diference between the take-o↵ and landing points can also be observed (cfr. Figure 5(b)). 

Before all flight lines were flown, a total of 31 control points (CP) which consisted of a CD glued at the center of a red cardboard of 31.5 cm x 50 cm (Figure 7a) were distributed throughout the bottom of the crater and on the lower region of the crater's rim-where the quadcopter was launched-as shown in Figure 7. The coordinates of each CP were determined with a Trimble R7 on Real Time Kinematic mode with a horizontal precision of 0.8 cm and 1.5 cm on the vertical and set to only register coordinates with a maximum Position Dilution of Precision (PDOP) of 4 (PDOP≤4). The time required to complete the field work, once on site, was of eight hours with a team of three people.

Figure 7 A total of 31 targets were used as ground control points in this study: (a) a red cardboard of 31.5 x 50 cm with a CD glued at its center was used for each target, (b) the targets were distributed throughout the domain with 26 targets located on the crater's bottom and 6 on the lowest part of the crater'stop. 

RESULTS AND DISCUSSION

Generation of Digital Surface Model

The camera was set to register the images in raw mode and transformed to 8-bit per channel tifs, with an approximate size of 37 MB per image. These tifs were processed with JetPhoto Studio (http://www.jetphotosoft.com), which matches the coordinates registered by the GPS watch with the acquisition time of each image, and writes the corresponding coordinate on each image's EXIF header. These coordinates are used by PhotoScan to speed up the image matching process, as the images are first aligned according to their coordinates and to set-up the coordinate system. With this alignment, Photoscan creates a cloud of matching points; if the matching points are correct, then a dense point cloud is generated either using the full resolution of the images (ultra high setting) or lower resolutions (high, medium and low), which downsample the original images by an order of four for each lower resolution. The final step consists of creating a three dimensional mesh based on the dense cloud, from which a Digital Surface Model and orthophoto can be exported as GeoTIFs; the dense cloud point and the 3D mesh can also be exported in different for mats for further processing in cloud point related software. For more details on the algorithms used in Photoscan, we suggest the review of ^{Verhoeven (2011)}.

Of all the acquired images, a total of 718 pictures were used in PhotoScan (Figure 8), because the pictures taken during launching and landing were discarded. Photoscan was used on a workstation running Debian 8.1 with two six-core Intel Xeon E5645 processors, 96 GB of RAM, and two GPU units (GeForce GTX 750 Ti and GeForce GTX 580). After the images were imported into Photoscan, they were calibrated through the use of Agisoft lens, which determines the camera's parameters based on several pictures of a chess-like pattern. The images were then aligned using the coordinates of the EXIF header, producing a sparse point cloud. This sparse point is used to verify the image alignment, after which a dense point cloud is generated. The dense cloud can be generated with four different settings: ultra high, high, medium, low and lowest. The ultra high setting uses the original images, while each lower setting downscales the original images by a factor of 4. On this work, the ultra high setting was used, requiring a total of nearly 29 hours and 56.2 GB of RAM, with PhotoScan's OpenCL option enabled. The ultra high quality generated dense point cloud comprised a total of 823,992,652 points, as shown on Figure 8 along the pictures used and the location of the CPs. The final step is to generate a 3D mesh of the model, for which PhotoScan uses either a Height or Arbitrary Height surface reconstruction; for topographic analysis and orthophoto generation, the Height Field Surface option is recommended in the Photoscan manual, and was the one used in this work. After the 3D mesh was created, both the DSM and the orthophoto were exported as GeoTIFs and imported into the GRASS GIS (^{GRASS Develoment Team, 2015}) in their native resolution, which was 4.7 cm. The resulting DSM and slopes of the crater's bottom are shown on Figure 9 as a composition with the generated orthophoto and shaded relief.

Figure 8 Aligned photographs on Photoscan Pro with the generated dense point cloud on the background. The GPS coordinates registered by the portable GPS attached to the quadcopter were embedded on each photograph's exif field in order to speed-up the alignment process. 

Figure 9 Derived products of the Digital Surface Model obtained in this work: (a) elevation, (b) Slope. Both maps are combined with the final orthophoto overlaid on a shaded relief of the final DEM (sun azimuth=225°, elevation=30°). Coordinates are in meters, UTM-14. 

Accuracy assessment

To quantify the errors of the generated Digital Surface Model at a resolution of 4 cm, we used the 31 Ground Control Points (GCPs), obtaining a horizontal RMSE=3.3 cm and a vertical RMSE=1.8 cm, with means of 2.6 cm and -0.3 cm, respectively; pixel-wise, the RMSE=0.3 pixel, which is an acceptable value (e.g. RMSE<0.5 pixel). The distribution of both horizontal and vertical errors is shown on Figure 10, along with their respective boxplot and histogram. Horizontally, nearly one third of the GCPs (10 out of 31) have an error between 1-2 cm (Figure 10a), with a maximum error of 7.5 cm. On the vertical, the residuals were determined as R = DSM − GPS, with residuals between 0-1 cm for 12 out of the 31 GCPs (Figure 10b), with one point located at the crater's rim registering a nearly -6.5 cm difference (with negative values representing lower elevations on the DSM than those measured with the RTK-GPS). This point also registered the largest horizontal error, with nearly 7.5 cm (Figure 10a); however, this difference is not observed on the other four points located on the crater's rim: horizontally, these four points had residuals below 3 cm; vertically, the two points located westward had residuals below 1 cm, while the two located eastward had residuals of -1.3 cm and -0.9 cm. The low accuracy on the crater's rim for that particular point was caused by a low photograph coverage on that area (Figure 8), as that is where the quadcopter turned to continue to the next flight line. However, it should be kept in mind that this point is located outside the domain of interest and that these five targets were placed to verify the accuracy of the generated model. On the crater's bottom, the largest horizontal error was of 7 cm (for two points), with three points showing errors of 6 cm; on the vertical, the largest error was of -4.5 cm; these errors were caused by the the high reflectance of the sediments.

Figure 10 Errors obtained on the ground control points: (a) horizontally and (b) vertically. The spatial distribution of errors is shown on top with their corresponding histogram and box plot. Negative values on the vertical indicate that the obtained DEM had lower elevations than those registered with the RTK GPS. The background image is the final orthophoto obtained in this work. Coordinates are in meters, UTM-14. 

To further assess the vertical accuracy of the elevation model obtained with Photoscan, we used the elevation measurements registered at 119 randomly located points, (Figure 11), with positive values indicating that the generated DSM had higher values than those measured with the RTK-GPS (i.e., residual=DSM-GPS). The vertical residuals of these random points have an RMSE=13.9 cm, with a mean=3.2 cm. As can be seen on the histogram of the residuals, 29% of the random points (35 out of 119) have residuals between 0-10 cm, with 19 points (16%) showing residuals between 0-5 cm. The distribution of the residuals (Figure 11) does not show a systematic error, but a random one, as both positive and negative residuals appear close to each other. To improve the accuracy assessment, an option would be to place targets similar to those used for the Ground Control Points (Figure 7a) but with different color (e.g., red for the GCPs and blue for the assessment points), to ease their location on the photographs.

Figure 11 Vertical errors observed on random points. Negative values on the vertical indicate that the obtained DEM had lower elevations than those registered with the RTK GPS. The background image is the final orthophoto overlaid on a shaded relief of the final DEM (sun azimuth=225°, elevation=30°). Coordinates are in meters, UTM-14. 

The detail of the generated 3D model is satisfactory, as it adequately represents the structures that have appeared due to the deformation of the lake's bottom, as clearly seen on Figure 12, where changes in fracture density, crestal grabens and mud injection domes can clearly be seen. With this methodology we were finally able to identify all the structures that have appeared in Parangueo and further work will be undertaken to determine the variation of its deformation rate.

Figure 12 Shaded relief (sun azimuth=225°, elevation=30°) of Parangueo, showing the fractures and domes that have been created due to its active deformation and mud injection. The shaded relief highlights the pronounced changes in fracture density along the fault scarp, the crestal grabens and the mud injection domes, with the large domes displaying both radial and concentric fracture/fault patterns. 

CONCLUSIONS

We presented a methodology to develop an ultra-high resolution Digital Surface Model (DSM) for a dry bed lake that shows a myriad of structures caused by active deformation (by both subsidence and uplift), such as escarpments, fractures and domes. Through the use of a small quadcopter, a point-and-shot camera, and a commercially available software that uses Structure from Motion algorithms we were able to develop both a DSM and an orthophoto at a resolution of 4.7 cm for a 1.44 km² area, with acceptable RMSE values (3.3 cm on the horizontal and 1.8 on the vertical). Rincón de Parangueo proved to be an excellent study site due to its complex topography and high reflection material to test our methodology, which produced a three dimensional model of unprecedented detail, showing all structures (fractures, domes and escarpments) that have been created due to the active deformation of the dry bed lake. The main advantages of the presented methodology are repeatability and reliability, which can allow the quantification of deformation rates in the study area through the comparison of DSMs acquired at different dates. This methodology is recommended from small to medium sized projects, where the costs of other type of data hinder their acquisition (i.e. LiDAR); although its application to large scale projects is feasible, the main restrictions of this methodology are: flight altitude of the UAV, restrictions on UAV operations, and the computing capacity available to process all required imagery.