(PDF) Sampling Designs for Selecting Points from Digital Surface Models to Create Digital Terrain Models

Sampling Designs for Selecting Points from Digital Surface Models to Create Digital Terrain Models Philip Fayad Apostolos Papakonstantinou Phaedon Kyriakidis Cyprus University of Cyprus University of Cyprus University of Technology Technology Technology Department of Civil Department of Civil Department of Civil Engineering and Geomatics Engineering and Geomatics Engineering and Geomatics School of Engineering and School of Engineering and School of Engineering and Technology Technology Technology 30 Archbishop Kyprianou Str., 2-8 Saripolou Str., Achilleos II 2-8 Saripolou Str., Achilleos I Limassol, Cyprus Bldg., 3nd floor, Bldg., 2nd floor,

Limassol, Cyprus Limassol, Cyprus

Abstract Recent advances in Unmanned Aerial Vehicle (UAV), in photogrammetric software and the miniaturisation of sensors lead topographers to embed the Unmanned Aerial System - Structure form Motion (UAS-SfM) pipeline to field work for earth measurements and DTM creation. This study attempts to develop a tool and examine a methodology to produce Digital Terrain Model (DTM) from a Digital Surface Model (DSM) created following UAS-SfM pipeline, while investigating the effect of known elevation point’s spatial distribution during interpolation. The proposed methodology aims to automate the DTM production process and minimize the production time, while maintaining a high-quality result and a low-cost approach. A geo-processing tool named “RS_Sampling Tool” is developed in order to extract known height points updating altitude information from the DSM using two sampling design schemas, random and stratified. The Inversed Distance Weighting (IDW) interpolation method was selected to create the DTM’s from six scenarios, in which the elevation points are differentiating in quantity and spatial distribution. The DSM used was created from a UAS flight realised on 25th of July 2017 to map Vrisa settlement on the island of Lesvos in Greece after a catastrophic earthquake that took place. The main conclusions of this study include: a) known altitude height points should be derived from bare ground areas of the DSM, b) the more points being used from the DSM the more accurate the final DTM is and c) a random sample distribution is much more likely to yield a DTM with low accuracy. However, a random based distribution can yield a high accuracy DTM, by including high-precision GCP’s. Keywords: Sampling Design; DTM; DSM; interpolation; methodology; elevation points. 1 Introduction interpolation methods, mainly because of its simplicity and efficiency. However, in order to obtain the best results a Regular and irregular grid, contouring and the sectional relatively dense elevation sampling distribution is required diagram are the common graphic representations (Gold, (Watson and Philip, 1985). 2005). Digital Terrain Models (DTM’s) are produced using Due to the very recent advances in the fields of computer elevation sampling points distributed in space, in a way that is vision and photogrammetry and in combination with the representative of the Area of Interest (AoI) and they can be improvements in data processing power, orthophoto maps and measured using various altitude recording devices like GPS, DSMs can be easily produced by terrestrial and/or aerial high- Lidar, etc. (Heywood et al., 2006). In this way it is possible to resolution 2D imagery. A plethora of scientists encourage the apply spatial interpolation methods and generate elevation use of Unmanned Aerial Vehicle (UAV), according to Adams values throughout the whole AoI, creating a DTM, a Digital (2011) and Papakonstantinou (2019), UAVs have been Elevation Model (DEM) or a Digital Surface Model (DSM) utilized with great potential following the 2009 L’Aquila, (Gold, 2005; Watson and Philip, 1985). In general, spatial 2010 Haiti, 2011 Japan and 2017 Vrisa earthquakes and each interpolation is the process of calculating unknown values in event presented different opportunities and lessons that will space, using some sampling observations (values) and it refers mould the promising future of UAV usage for imagery to continuous phenomena throughout the space (Heywood et collection in disaster management and monitoring. al., 2006). The terms “DEM” and “DTM” are more specific Furthermore, an extensive study was implemented to monitor and refer only to digital altitude representation of the earth’s and to map Vrisa village damage assessment as a post- surface without any additional elements, while DSM’s refer to earthquake process (Soulakellis et al., 2018). surfaces that also contains other physical or anthropogenic Recent technological advantages make UAV-based elements such as structures, vegetation, etc. photogrammetry highly suitable for surveys in a geo-hazard In this study the Inversed Distance Weighting (IDW) context, as in a post-earthquake scenario, and its advantages technique is used, one of the most popular deterministic may be summarized as follows: a) safety: no risk for AGILE 2019 – Limassol, June 17-20, 2019 operators; b) possibility to survey inaccessible zones; c) high- resolution photographs; d) speed of survey and elaboration; Figure 1: Methodological steps. and e). repeatability and economic convenience (Dominici et DSM al., 2017; Kavroudakis et al., 2018). production This study aims to develop a semi-automated method for selecting points from a DSM in order to produce a very Object-Based detailed in terms of resolution and earth’s surface ground Classification variation DTM. A UAV flew over the Area of Interest (AoI) and the DSM is produced following the Unmanned Aerial Sampling Design System - Structure form Motion (UAS-SfM) pipeline. A total of six sampling scenarios were implemented to derive ground DTM heights from the UAS-DSM. The heights derived from the production DSM were used for the creation of new a DTM, one for each scenario. The DTMs produced are compared in order to On the 25th of July a UAS data acquisition campaign from investigate the effect of known elevation point’s spatial the University of the Aegean took place for mapping Vrisa distribution during interpolation. More specific in the DTM Settlement at village spatial scale. In this campaign the comparison the spatial distribution of the elevation points is following spatial data were collected: a) 148 Ground Control being investigated by varying the quantity and the randomness Points (GCP’s) using a Real Time Kinematic (RTK) system of the elevation points used in the spatial interpolation and b) 1000 very high-resolution aerial images. The UAS flew process. In the following chapters the proposed methodology at low altitude (65 m), capturing high resolution images with is analysed. 80% overlap and 80% sidelap. Following the UAS-SfM pipeline, the orthophoto map (fig.2) and the DSM (fig.3) created having a spatial resolution of 3 cm and 5 cm, 2 Materials and Methods respectively. Traditionally the creation of DTM for an area is implemented Figure 2: The Orthophoto map produced. using spatial interpolation to ground measurements. In this process the calculation of unknown values based on some sampling observations applies to DTM creation. 2.1 Area of Interest On 12 June 2017 (UTC 12:28:38.26) a magnitude Mw 6.3 earthquake occurred offshore Lesvos Island in SE Aegean Sea, which was widely felt to the citizens of the island. Most substantial damage was reported to the village of Vrisa located to the south-eastern coast of the island. The study area covers all the damaged part of the village, approximately 0.3Km2. Figure 3: The Digital Surface Model produced. 2.2 Methodology The workflow of the semi-automated methodology presented in this study (fig.1) consists the following steps: 1. DSM production: a DSM of the AoI is produced following UAS-SfM pipeline, 2. Object-based Classification for the automatic identification and selection of pure ground polygons on the UAS-DSM, 3. Selection of elevation points within the polygons selected from previous step according to the parameters (RS Sampling geoprocessing Tool creation), 4. Spatial interpolation using the IDW approach for the DTM production. From the constructed DSM it is crucial to clearly identify bare ground areas in order to select and extract elevation values from these areas. This can be achieved using classification methods. Two are the most important: a) the pixel-based, in which classification is done according to the ground spectrum differences and b) the object-based method, AGILE 2019 – Limassol, June 17-20, 2019 in which classification is relying not only on the spectral scenarios. That was achieved through the Sampling Design characteristics of ground, but also consider geometric and process and using the random or the stratified sampling structural information. Many studies, like the one method. The number of points used was not the same in each implemented by Chen (2009), had proven the efficiency and of the six scenarios thus, consisting: a) a standard number of the advantages of the Object-based classification method over 626 points that represent the centroid of the polygons of small the traditional pixel-based method especially in reducing areas (areas < 10m2) and b) in some scenarios, 55 (out of the errors/noises and “salt and pepper” phenomena. Figure 4 total 148) additional high precision (RTK) GCP’s. The display the process followed in order to define the elevation distribution of the 55 GPS points followed a pattern that cover points in each one of the six scenario. all the study area. The 35 were selected visually according to To define polygons of the bare ground areas in the DSM, the terrain variations of the study area and the remaining 20 the Object-based classification approach was implemented, according to their height value (where the average height using the Feature Extraction tool in ENVI 5.0 (Feature value of the area is exceeded). Extraction Module User’s Guide, 2008). This led to 1250 In this study, six point sampling scenarios were applied polygons of various sizes containing a total of 25+ million (table 1) using the random and stratified sampling methods. In elevation points. The 626 of the total 1250 polygons were D, E scenarios the preselected 55 additional GCP’s were polygons with area less than 10 m2 (labeled as small areas). added to the DTM creation. The sampling parameters used The result of this classification process was polygons that were the following: a) number of points in each polygon, b) contained a large number of elevation points. Following the minimum allowed distance between points and c) Cell size. methodology (fig.1), in order to select the appropriate amount The number of points distributed in each polygon was set to of elevation points needed and reduce the total number so that 100 at scenario A while at scenarios B and D was depended those can be used in the spatial interpolation process, the next from the point height range. The minimum allowed distance step was to create the different elevation points selection between points was set as a fixed value (2m) at scenarios A, B Table 1: Elevation points selection scenarios. Sampling Sampling No of Points Interpolation Interpolation Scenario method parameters used method parameters Number of points in each polygon: 3,209 p = 2, fixed = 100, + 626 (of small A Random IDW Search radius Variable, Minimum allowed areas) Number of points = 16 distance between = 3,835 points: fixed = 2m Number of points in each polygon: Variable = 1,557 p = 2, depending on + 626 (of small B Random IDW Search radius Variable, height range, areas) Number of points = 16 Minimum allowed =2,183 distance between points: fixed = 2m 7,109 p = 2, Cell height = 3m + 626 (of small C Stratified IDW Search radius Variable, Cell Width = 3m areas) Number of points = 16 = 7,735 Number of points in each polygon: 1,612 Variable = + 55 (GPS) p = 2, depending on D Random + 626 (of small IDW Search radius Variable, height range, areas) Number of points = 16 Minimum allowed = 2,293 distance between points: fixed = 2m 7,109 + 55 (GPS) p = 2, Cell height = 3m E Stratified + 626 (of small IDW Search radius Variable, Cell Width = 3m areas) Number of points = 16 = 7,790 25,569,387 p = 2, All elevation + 626 (of small Z - IDW Search radius Variable, points are used areas) Number of points = 16 = 25,570,013 AGILE 2019 – Limassol, June 17-20, 2019 and D. At scenarios C and E the parameter cell size was set to The “RS Sampling Tool” developed implements the Sampling 3x3m in height and width respectively. Finally, at scenario Z design calculations and exporting all the resulting elevation all elevation points delivered from the DSM were used. In all points, according to the parameters set. The last step in the scenarios, the IDW method was used, with the “Search proposed methodology is the production of the DTM using radius” and “Number of points” as predefined variables the spatial interpolation process. having values 2 and 16 respectively. Figure 4: The elevation points selection process (Sampling Design) for the interpolation and the DTM production. 3 Results Figure 5: The resulting DTM surface in each scenario. The resulting surfaces produced from the proposed methodology were of small differences (fig.5). All DTMs produced, were having heights variation of 25.503m (scenario C and E) to 73.373m (scenario B). Surface Z (fig.6) was an exception having the lowest height (25.390m) as well as the highest (74.019m) height values among all surfaces. Comparing the six elevation surfaces, scenario Z has resulted to a much rougher surface and that is due to the fact that the total of the 25 million elevation points was used. The height variation in all the scenarios implemented is depicted in Table 2. Furthermore, a statistical comparison of the DTMs was implemented using the min-max value matrix (table 2), the RMSE table (table 3 and fig.7) and the correlation matrix (table 4). To investigate and visualize the differences between the DTM of scenario Z with all the other scenarios, we used the cut fill tool in ArcGIS (fig.8). Table 2: Min and Max height value in each scenario. Scenario MIN (m) MAX (m) A 25.629 73.468 B 25.919 73.373 C 25.503 73.609 D 25.740 73.374 E 25.503 73.609 Z 25.390 74.019 As for the accuracy of the produced DTMs, table 3 depicts the total RMSE errors, which are calculated using the remaining GCP’s in each scenario, out of the total 148 GCP’s. AGILE 2019 – Limassol, June 17-20, 2019 Figure 6: The resulting DTM of scenario Z, where all 25 million elevation points are being used. Table 3: The calculated RMSE in each scenario. Scenario RMSE RMSE Number of Elev. (all 148 GCP’s used) (148 -55 = 93 GCP’s used) points A 1.007m ----- 3,835 B 1.108m ----- 2,183 C 0.620m ----- 7,735 D ----- 0.897m 2,293 E ----- 0.590m 7,790 Z 0.377m ----- 25,570,013 Figure 7: The number of elevation points used in the DTM generation in correlation with the RMSE that resulted in each scenario. 1.20 1.108 1.007 1.00 0.897 0.80 RMSE (m) 0.620 0.590 0.60 0.377 0.40 0.20 0.00 2,183 2,293 3,835 7,735 7,790 25,570,013 Number of elevation points Table 4: Correlation Matrix where values range between 0 and 1. Zero value means absolutely no correlation between the scenarios while a value of one indicate exact same surfaces. Scenario A B C D E Z A 1 0.99974 0.99981 0.99975 0.99981 0.99907 B 0.99974 1 0.99951 0.99999 0.9995 0.9988 C 0.99981 0.99951 1 0.99954 1 0.9993 D 0.99975 0.99999 0.99954 1 0.99954 0.99883 E 0.99981 0.9995 1 0.99954 1 0.9993 Z 0.99907 0.9988 0.9993 0.99883 0.9993 1 AGILE 2019 – Limassol, June 17-20, 2019 Figure 8: Volume changes (differences) between all the DTM surfaces (scenarios A to E) from the DTM of scenario Z. Red color indicate Gain (addition), grey represent unchanged areas and blue indicate Loss. 4 Conclusions Finally, as the quantity of elevation points that inserted to the DTM generation increases the greater the resulted In this study, a semi-automated methodology for creating a accuracy of the final product (fig.7) DTM using Z values from points delivered from a DSM is proposed and tested on a real AoI. The implementation of six point sampling scenarios was realized using as a ground truth References values 148 high precision (RTK) GCP’s. Thus, the effect of point sampling quantity and point distribution on the quality Adams, S. & Friedland, C. (2011) A survey of unmanned of the final DTMs was examined. aerial vehicle (UAV) usage for imagery collection in disaster Comparing the total accuracy of the DSM from the UAS- research and management. 9th International Workshop on SfM pipeline which is 0.645m, with the accuracy of the DTM Remote Sensing for Disaster Response, 15-16 September produced in each scenario (table 3), it is clear viewed that 2011. California, USA. scenarios C, E and Z are creating DTMs of higher quality. Must be noted that scenario Z has nearly as twice as accuracy Chen, J., Xiang, L. & Jingjing, L. (2009) The Study of Object- in comparison to the other two scenarios. Additionally, a close Oriented Classification Method of Remote Sensing Image. look at scenarios A and B, RMSE results leads us to the First International Conference on Information Science and following conclusion. Α random sample distribution is much Engineering, ICISE 2009, 26-28 Dec 2009. Nanjing, China. more likely to generate a DTM with low accuracy rather than pp. 1495–1498. a one with high, mainly because the random sampling. Moreover, the ground altitude variations are not considered as Dominici, D., Alicandro, M. & Massimi, V. (2017) UAV a parameter during the selection of the elevation points. photogrammetry in the post-earthquake scenario: case studies However, the results of scenario D indicate that a random in L'Aquila. Geomatics, Natural Hazards and Risk. [Online] based distribution can generate a high accuracy DTM as long 8(1), 87-103, Available from: doi: as a number of high-precision GCP’s is included. 10.1080/19475705.2016.1176605 [Accessed 21th Jan 2008]. Additionally, scenarios C and E tend to produce DTM having similar surfaces (see table 4) and maximum correlation value ENVI image analysis software. (2008) ENVI Feature of 1. This leads to the conclusion that the additional 55 high Extraction Module User’s Guide. [Online] Available from: precision (RTK) GCP’s had no impact when used to the http://www.harrisgeospatial.com/portals/0/pdfs/envi/feature_e stratified sampling method. xtraction_module.pdf [Accessed 20th February 2018]. AGILE 2019 – Limassol, June 17-20, 2019 Heywood, I., Cornelius, S., & Carver, S. (2006). An introduction to geographical information systems. Harlow, Pearson Prentice Hall. Kavroudakis, D., et al. (2018) Efficiency and Effectiveness Approaches in Spatial Data Collection of Vrisa after Lesvos Earthquake. International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Gi4DM, 18–21 March 2018, Istanbul, Turkey. 42(3W4), 277–281. Li, Z., Zhu, C. & Gold, C. (2005). Digital Terrain Modeling: Principles and Methodology. Washington, D.C., CRC Press. Soulakellis, N., Tataris G., Papadopoulou E., Chatzistamatis S., Vasilakos C., Roussou O., Papakonstantinou A., Kavroudakis D., Kalloniatis C. (2019) Synergistic Exploitation of Geoinformation Methods for Post-Earthquake 3D Mapping and Damage Assessment, In: Orhan, A., Chandra, M., Sunar, F., Tanzi, T.J. (eds.) Intelligent Systems for Crisis Management: Gi4DM 2018 (Lecture Notes in Geoinformation and Cartography). Springer. Soulakellis, N., et al. (2018) Synergistic Exploitation of Geoinformation Methods for Post-Earthquake 3D Mapping of Vrisa Traditional Settlement, Lesvos Island, Greece. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Gi4DM, 18–21 March 2018, Istanbul, Turkey. 42(3W4), 491-498. Watson, D.F. & Philip, G.M. (1985) A Refinement of Inverse Distance Weighted Interpolation. Geo-Processing, 2(4), 315– 327.