Quality Assessment of Globally Available Digital Elevation Model (DEM) for Mountainous Region

1. Introduction

The topographic representation in the form of Triangulated Irregular Network (TIN), spot height (x, y, z) and contour constituting  the Earth’s surface is called Digital Elevation Model (DEM) (Forkuor & Maathuis, 2012). DEMs are considered as preconditioning datasets in geospatial applications including extraction of terrain features for geomorphological applications, water flow direction for hydrological and hydrodynamic modeling, landscape dynamic and ecosystem modeling (Dong et al., 2015; Ioannidis et al., 2014). With the availability of DEMs from SRTM, ASTER GDEM, and ALOS PALSAR platform, the inventorization of topographical parameters and its derivatives eased the uses of robust and expansive ground-based spot height survey and UAV based InSAR using ground control point (Sze et al., 2015). SRTM DEM operates in RADAR C-band microwave electromagnetic spectrum, launched in 2000 by NASA Earth observing system. Beside this, it also uses two sets of sensor vertically downward and off nadir directions to generate the stereo-pair imagery for the construction of DEM through InSAR process (Rabah et al., 2017). In addition, ASTER, GDEM is an optical dataset which provides global coverage with 30m spatial resolution on Terra Spacecraft launched in 1999 by  NASA’s Earth Observing system in joint collaboration with Japan’s Ministry of Economic Trade and Industry (Rabah et al., 2017; Yap et al., 2018). Similarly, the Terra spacecraft was designed with off nadir and backward viewing sensor onboard providing dual band off-nadir and backward viewing band producing stereo-pair imagery which is used for DEM generations covering 83ºN and 83ºS (Dong et al., 2015). ASTER GDEM is georeferenced to WGS 84 and EGM96 model with the geographic coordinates in latitude and longitude from the “sensor in location tagged image format” (GeoTIFF) (Chirico et al., 2013). ALOS PALSAR satellite was equipped with the along-track 2.5m resolution PAN sensor for Stereo mapping (PRISM), 10m resolution Advanced Visible and Near-Infrared Radiometer (AVNIR) and L-band Polarimetric Phased Array Synthetic Aperture Radar (PALSAR) launched by JAXA’s Tanegashima Space Center in southern Japan in 2006 orbiting sun synchronously (Bignone and Umakawa, 2008).

Perhaps, the radar beam has the capability to penetrate  cloud and vegetation canopy through certain extent, however in sparse vegetation, the radiation gets scattered back from the surface to the sensor  leading to distortions such as layover, foreshortening and shadow (Yap et al., 2018). Whereas under the dense thick vegetation, the c-band electromagnetic radiation and optical based ASTER GDEM doesn’t penetrate the earth surface, and therefore doesn’t reflect the actual backscatter signal of the earth surface (Gou-an et al., 2001). This exerts a necessity for geospatial analysts, researchers and planner to examine the qualities of DEMs prior to the development of applications. Numerous authors have assessed the qualities of DEMs in the past. Chirico et al. (2013), has evaluated the accuracy of ASTER GDEM v1 and v2 using geomorphic derivatives of DEMs; land use landcover, slope, relief and topographical terrain. The research found that the ASTER GDEM v2 provides a better accuracy as compared to ASTER GDEM v1 in heterotypic land use land cover features, highland and bare land. Blanchard and Rogan (2014) has assessed the geomorphic change analysis based on ASTER and SRTM DEM’s and observed that ASTER GDEM provides higher accuracy over densely populated forest area and lower over agricultural and grasslands. However, SRTM provides high accuracy over bare land and built-up areas and lower over dense vegetated areas. ALOS PALSAR DEM have been found to be accurate in term of horizontal offset for hydrological applications especially for the derivation of drainage line with coefficient line of correspondence (CLC) of 0.91 demonstrating the strong correlation between the visually interpreted drainage line (Niipele and Chen, 2019). The reduction in topographical effects has increase the vertical accuracy of 0.98 using assisted offset tracking procedure of Muztagh Ata mountain (Yan et al., 2013) and the similar methods was applied to enhanced the vertical accuracy when topographic correction has been made by interferometric SAR process using stereo-pair imagery from ALOS PALSAR instrument, however degradation, mis-registration and improper estimation of backscatter signals was vividly observed while topographic correction was made with coarser spatial resolution of external global DEMs (i.e. ASTER and SRTM) (Das et al., 2014). Most of the previous studies focused on vertical error using ground control point (herein refer as GCP) (Baldi et al., 2010; Dong et al., 2015; Yap et al., 2018), existing toposheet (Saran et al., 2010; Sharma et al., 2010) and high-resolution DEM generated from ortho image using stereoscopic process (Baldi et al., 2010),  however, the DEMs  are accountable for both vertical and horizontal errors (Rabah et al., 2017). Therefore, the study emphasized on evaluation of both vertical errors using ground control points (GCP) and horizontal displacement using Topo-DEM for the globally available DEMs.

Previously, numerous image matching algorithm have been applied for various applications across the globe (Heid and Kääb, 2012; Zhao et al., 2006). The image matching algorithm has been discussed and demonstrated with the level of accuracy in context to glacier flow and ground displacement by (Heid and Kääb, 2012). The study evaluated the accuracies of six different type of image matching algorithm; Normalized Cross Correlation (NCC) techniques, Cross Correlation using Fast Fourier transform (CCF), Phase Correlation using FFT (PC), Cross Correlation using FFT on orientation Image (CCF-O), phase Correlation on FFT on orientation image (PC-O) and Cosi-correlation techniques. Among those, for small search window with narrow and heterogeneous landform, NCC is found to be higher accurate while CCF-O perform best in the context of poor and good visual context in larger search window size. Therefore, as a pre-requisite to geospatial applications, the assessment on quality of DEMs is significant and necessity in the mountainous region prior to implementation. Therefore, the study aimed at approaching NCC, RMSE statistical model and sub-pixel offset methods to assess the qualities of DEMs in the mountainous region.

2. Materials And Methods

The freely available global DEMs i.e., SRTM, ASTER GDEM and ALOS PALSAR were used in the study. 5m contour topography surveys were conducted by Phuentsholing City Corporation through conventional ground based total station (TS-06) in 2016 for the entire city of Phuentsholing was collected to derive the topographical DEM (herein refer as Topo-DEM) which will be reference to evaluate the other DEM’s. SRTM v3 and ASTER DEM v2 of 30m spatial resolution was downloaded from USGS earth explorer and 12.5m ALOS PALSAR DEM from Alaska Vertex website. All the DEM’s were converted from UTM 46N projected coordinate system with respect to WGS 84 datum to DRUKREF 03 Bhutan National Grid with respect to D-Bhutan National Geodetic Datum coordinate system and analyzed based on local coordinate system. In order to determine the vertical accuracies of the DEM’s 100 GCPs using GPS were collected around the study area. Several researchers have evaluated the qualities of DEMs with vertical RMSE statistic across the globe demonstrated in Table 2. The sources of various material used in the study and typical representation of different DEMs are shown in Figure 2 and Table 1.

2.1 Study Area

The study area consists of undulating topography with structural geological features which fall under sub-Himalayan region in Bhutan. Phuentsholing is geographically bounded by 26°51'41.51''N and 89°23' 0.95'E with an elevation of 182m above mean seal level. The region is bisected by the parallel drainage pattern of Om Chhu River and valley formation by Ammo Chhu river toward southern part of the study area. The region was covered with rugged terrain features, sparse vegetation cover, and abrupt rise in topography from Indian flood plain along the Ammo Chhu River. Purposely, Phuentsholing area was selected for the analysis of DEM qualities, since the area consist of both plain and undulating topography to compare RMSE of both horizontal and vertical offset with respect to Topo-DEM. The study area map with elevation ranges along with randomly distributed GCPs is shown in figure 1.

2.2 Methods

2.2.1 DEM Generation from Topography maps (Topo-DEM)

The topographic survey was conducted in 2016 by Phuentsholing City Corporation. 5m contour survey data from topography map was used for generation of digital elevation model of 5m spatial resolution. However, the spatial resolution of SRTM and ASTER GDEM is 30m. Therefore, derived Topo-DEM and 12.5m ALOS PALSAR DEM was further converted to 30m spatial resolution using bilinear convolution techniques to rationalize with other DEMs. 100 randomly distributed samples were gathered using kinematic hand held GNSS receiver to compute the vertical RMSE error of each DEM with respect to Topo-DEM.

2.2.2 Root Mean Square Error (RMSE)

The RMSE, in geospatial domain with regard to the present research defined as difference in elevations of the global DEMs with respect to the high-resolution ground base Topo-DEM derived from topographic map which is represented by equation 1.

\(RMSE = \sqrt{\sum{ {(P_i-Q_i)^2} \over n}}\)                                  (1)

where, ‘P_i’ refers to the elevation values for Topo-DEM, ‘Q_i’ denote the elevation of other tested DEMs and ‘n’ represent the number of randomly distributed samples in the study area. The RMSE provide the elevation difference between the tested and high-resolution corrected Topo-DEM in term of vertical accuracy, while the normalized cross-correlation (NCC) between different DEMs with respect to the Topo-DEM, determines the correlation among the DEMs and sub-pixel based approach which evaluates the horizontal shift and displacement among the DEMs (Chirico et al., 2013; Elkhrachy, 2018).

2.2.3 Normalized Cross Correlation (NCC) and Sub-pixel Offset techniques

Cross correlation techniques have been used widely with the advancement of digital transformation from analogue system. NCC is the tool that measure the degree of similarity between the two-digital image using correlation coefficient algorithm (Debella-Gilo and Kääb, 2011; Rao et al., 2014). Besides, cross correlation, NCC also evaluates the degree of similarity or dissimilarity between the two images. The normalized correlation can be represented by the formula given in Equation 2.

\(NCC' {_(x,y_)} (t) = {NCC{_(x,y_)(t)} \over {\sqrt NCC_(x,x_) (0) NCC _(y,y_) (0)}}\)

                    (2)

where \(NCC' {_(x,y_)} (t)\), represent the normalized correlation at acquisition time ‘t’, \(NCC' {_(x,y_)} (t)\)  denote the correlation of one image at location x, y at image acquisition time ‘t’ and \(NCC' {_(x,x_)} (0)\) .\(NCC' {_(y,y_)} (0)\)on another image at time ‘0’. The \(NCC' {_(x,y_)} (t)\)  values ranges from -1 to +1 where, +1 represent the both the image are exactly correlated, 0 represent the image are uncorrelated and -1 represent similar to +1 which is differ by negative sign only (Rao et al., 2014).

Debella-Gilo and Kääb (2011) and Heid and Kääb (2012) has assessed the velocity dynamic of the glacier flows and ground displacement with an assist from satellite imagery acquired at different time using different image matching algorithm namely; Normalized cross correlation techniques, Cross and Phase correlation on frequency domain using Fast Fourier transform on oriented and non-oriented image and COSI Correlation techniques. The research demonstrates that the area with heterotypic geographical features provide greater accuracies by NCC technique while in the case of homogenous geographical surfaces, is of less accurate compared to other unnormalized matching techniques. The general approach of the entire research is depicted in figure 3.

While, considering the horizontal error, the actual shift in pixel reference to the actual position of the pixel on the real ground is crucial and vital to be considered (Jin et al., 2014). The research also found out that the pixel offset in the DEMs are encountered in most of the global DEMs due to the process in georeferencing and DEM generation through photogrammetric technique. The sub-pixel offset measurement was undertaken with the help of Raster to Topo conversion tool in ArcGIS. The pixel of each resolution of \(30 × 30m\)  in vector format is extracted to determine the horizontal and vertical pixel offset of each DEMs with respect the Topo-DEM.

3. Results And Discussions

3.1 Descriptive Statistics

The descriptive statistics on the relative dispersion from the central of tendency (µ) was examined for the respective DEM using Coefficient of variation (CV) and standard deviations (SD). In general, Topo-DEM varied with an elevation from 185.21m to 593.44m with an average of 332.24m, 156.2m to 593.29m with an average of 320.71 for ASTER GDEM v2, 183.64m to 592.34m with an average of 336.34m for SRTM v3 and 139.31m to 557.39m with an average of 292.42m for ALOS PALSAR DEM, respectively. The unit less CV varied with the minimum of 0.29 for Topo-DEM and SRTM v3 followed by ASTER GDEM v2 (0.31) and ALOS PALSAR (0.33).

3.2 Computation of Root Mean Square Error (RMSE)

More than 100 randomly distributed ground truth samples over the study area were collected with the help of kinematic hand held GNSS receiver to assess the vertical erroneous of the globally available DEMs. The individual elevation of the DEMs was extracted to GCPs using “Add surface information values” tools in ArcGIS software. The elevations values were further exported to excel sheet for evaluation of RMSE of the individual DEMs. The RMSE of the DEMs are represented in table 2.

From the analysis, the RMSE of SRTM with respect to Topo-DEM is observed less than ASTER GDEM followed by ALOS PALSAR with 7.98m, 17.4m and 40.3m, respectively. Similar observations have been made as listed in Table 3 for undulated topographical region. The variation in RMS error at different region depends on the topographical features, land use types and geomorphological characteristics within the research area.

3.3 Computation of NCC

For NCC, the respective DEM files were converted to contour file with 30m interval. NCC was computed for both plain and hilly region to carryout comparative analysis with the reference Topo-DEM. With regards to the horizontal cross correlation, SRTM DEM follows the reference Topo-DEM, however with little shift in both ‘x’ and ‘y’ directions. The correlation coefficient between the reference Topo-DEM and SRTM DEM shows strong positive correlation with 0.995 as shown in figure 6.

Additionally, in context to ASTER GDEM, the contour doesn’t follow the reference DEM and mis-registration occurs in between with the horizontal (x) and vertical (y) shift. Further, mis-registration in contour can be observed in every interval of the ASTER GDEM with reference Topo-DEM. However, the strong correlation coefficient was also observed between the two DEM as shown in figure 6. From the observation made in figure 6, the NCC of DEMs demonstrates strong agreement in plain region as compared to hilly terrains.

3.4 Pixel Offset and Displacement Measurement

The pixel offset measurement was undertaken manually by converting the raster pixel into polygon using ArcGIS for all the DEMs reference to the Topo-DEM.  The measurement was taken along East-West and North-South directions. Besides, the displacement was also measured by the pixels displaced from the original reference pixel (Pixel value of Topo-DEM) to the pixel values of tested DEMs. The minimal offset of 0.48m toward east and 7.45m toward south was observed in the case of SRTM DEM while no displacement of pixel has been observed. Similar observation has been made for SRTM v3 by Jin et al. (2014) with half a pixel offset from the reference DEM. Additionally, for ASTER GDEM v2, the offset of 8.52m toward East and 9.06m toward North direction was observed with no displacement in plain region. In addition, for ALOS PALSAR DEM, the negligible offset of 1.5m toward south and 0.95m toward east was observed with 5 pixels displaced toward East and 10 pixels toward North with reference to Topo-DEM at plain area. The details offset and displacements are shown in figure 8 and table 5.

In the case of sloppy terrain, SRTM v3 follows the reference Topo-DEMs with negligible offset, while, both ASTER v2 and ALOS PALSAR observed the horizontal misregistration of an isoline of 30m from the reference Topo-DEM (Figure 5).

4. Conclusions

The globally available DEMs were converted to the local coordinate system DRUKREF03 National Grid. The kinematic hand held GNSS receiver was used for gathering ground control points (GCP’s). The normalized cross correlation was also computed to compare among the DEMs and observed with high correlation coefficient for all the DEMs value higher than 0.95. However, the 30m SRTM DEM v3 is better in such kind of area covered with relief, sparsely vegetated coverage and mixed with flood plain geomorphic features. The mis-registration are the common limitation observed in all the DEMs, which may be due to the process involved in interpolation, conversion and data manipulation in creating DEMs and its derivations.

The RMSE, horizontal and vertical shift and displacement of the DEMs with respect to the Topo-DEM has been assessed. From the analysis, it was observed that SRTM DEM v3 provide high accuracy than other DEMs with RMSE values less than 10m and horizontal shift of 0.5m toward East and 7.45m toward South direction. In addition, the RMSE of ASTER GDEM v2 was observed with 17.4m with horizontal shift of 8.52m toward East and 9.06m vertically toward North direction. However, ALOS PALSAR DEM was observed with RMSE 40.30m with negligible shift in horizontal and vertical direction but the pixel displacement of 5 pixels toward East and 10 pixels toward North from the reference Topo DEM.

Further, in continuation to the current study which is based on physical calculation on offset measurement and pixel displacement, an automatic code or a tool will be developed for image matching and derived the pixel offset, displacement and its correlation coefficient among the DEMs to reduce the RMSE and pixel offset to generate high accurate DEM to be applied on different purposes.

5. Funding Agency

This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.