Spatial Enhancement of DEM Using Interpolation Methods: A Case Study of Kuwait’s Coastal Zones

Digital elevation models (DEMs) are essential tools utilized in several branches of science, including environmental, geological, and geospatial studies. Unfortunately, high-accuracy DEM data such as LiDAR are not publicly available, and the coverage is limited. Therefore, the use of alternative methods, such as interpolation techniques (i.e., kriging, inverse distance weighting, radial basis functions), is greatly advantageous for the production of enhanced DEMs. The results of this study show that interpolated DEMs had minimal errors (RMSE = 1.44) with an increase of about 28% from the original DEM. However, the spatial resolution of interpolated DEM data was enhanced significantly by 83%. The deterministic interpolation methods provided more accurate estimations for producing DEMs in the coastal zones of Kuwait than geostatistical interpolation methods. The reference elevation data were collected using GPS and accurate topographic maps (1:25,000), and elevation points from the interpolated DEM were matched significantly (R 2 = 0.88; R 2 = 94, respectively). Given the lack of accurate DEM data, the interpolated DEM produced in this study are held in high regard and highly recommended for use in the coastal zone of Kuwait.


Introduction
Digital elevation models (DEMs) are a virtual representation of elevations where the terrain is established by three-dimensional coordinates [1]. Recently, DEM has become an essential tool in several branches of science, including environmental, geological, and geospatial studies. It is generated mostly in raster format using different sources, such as field measurements collected by GPS, topographic maps, and remotely sensed data [2,3].
The DEMs produced from these sources vary in terms of price, accuracy, availability, and sampling density [4]. In recent decades, the primary source of DEMs used in many types of research has been remote sensing due to its low costs, coverage, availability, and consistency with required spatial and spectral resolutions [5]. The vertical accuracy of remotely sensed DEMs varies depending on their source. LiDAR data, for instance, have a vertical accuracy of approximately 15 cm [6], whereas the vertical accuracy of DEMs from SRTM reaches ±16 m [7].
Numerous sources of error affect the quality of DEMs, but most factors can be attributed to data acquisition [8]. Unfortunately, accurate sources of DEMs, such as LiDAR, are not available in many places, especially in developing countries. They have limited spatial coverage, and require experts to process the data. Therefore, using alternative methods, such as spatial interpolation techniques, is necessary to produce DEMs with acceptable accuracy and resolution.
Spatial interpolation techniques are used to predict unknown values, such as elevations and air temperatures, when other values are known. The accuracy of the predictions is greatly increased as the number of known values increases [9]. The spatial interpolation could be conducted using deterministic or geostatistical methods. Deterministic interpolation techniques, such as inverse distance weighting (IDW) and radial basis function (RBF), create surface grids from the measured points based on similarity or degree of smoothing. Geostatistical techniques, such as kriging methods, utilize the statistical properties of measured points to generate a surface map [10].
This research aims to produce high spatial resolution DEM using the interpolation techniques to overcome the lack of highly accurate DEM data in the coastal zone of Kuwait and verify the quality of interpolated DEMs using ground controls points (GCPs). Several datasets were employed to produce high spatial resolution DEMs, including remotely sensed data, contour maps, and elevation points obtained from Global Positioning System (GPS). The DEM generated in this research could provide services in various environmental fields, particularly in climate change and associated issues, such as sea level rise.

Methodology
The elevation data used to generate a high spatial resolution DEM layer were obtained from DEM produced by the study [11] and contour maps with different scales. These spatial datasets underwent pre-processing as they were gathered from different sources with variations in accuracy and scale using ArcGIS application. The study [11] reduced the vertical accuracy of the ASTER Global Digital Elevation Model (GDEM) covering Kuwait from 17 m to 5 m. The elevation data from this DEM was employed in the present research.
Moreover, there were three topographic maps used in this research: one was a paper copy, so it was digitized, georeferenced, and registered, like other spatial datasets, to Universal Transverse Mercator, coordinate zone 39, and WGS 84. Elevation points obtained from two topographic maps with scales of 1:400,000 and 1:100,000 were utilized in the interpolation process, while elevation points from a higher resolution topographic map (1:25,000) were used to assess the vertical accuracy of the interpolated DEM.
Next, elevations datasets were subset to correspond to the coastal zones of Kuwait and then converted to points. The elevation points were interpolated to refine the spatial resolution to 5 m. After that, the DEMs produced by different kinds of interpolation techniques were evaluated using a cross validation method and correlation coefficient (r). The former method was performed by removing one elevation point and predicting this point with the model constructed using the remaining elevation points. The purpose of using cross validation was to evaluate the accuracy of the models in estimating unknown elevation values. For a best-fit model, the root mean square error (RMSE) should be as small as possible [12]. The latter method was used to calculate the association between the estimated elevation points and measured values of elevation points using correlation coefficient (r). Thereafter, the most acceptable interpolated DEM with lower RSME and higher (r) in each zone was selected. Then the interpolated DEM in each coastal zone were combined in one raster file. Then, the vertical accuracy of the interpolated DEM was examined by matching reference data (which were GCP obtained from both GPS and high-resolution topographic maps (1:25,000)), with corresponding elevation points in exact locations using linear regression coefficient (R 2 ), following the procedure outlined by [11]. These GCP data were selected clear of buildings to fit spatially with data of the interpolated DEM. Figure 1 illustrates the flowchart of the methodology.

Study Area
The State of Kuwait is located in the northwestern part of the Arabian Gulf. Its coastline extends 500 km, including its islands. Two of them (Boubyan and Failakah) occupy more than 95% of the total area of the islands ( Figure 2). Kuwait Bay is one of the most distinct features in the Kuwait marine environment. It is a semi-enclosed shallow body of water extending about 35 km inland. The coastal zone of Kuwait is a low-lying area, which makes it a prime location for the development of industries, residences, recreation, and tourism. The coast of Kuwait can be divided into three geographical regions: the north, which is non-urbanized; the south, including Kuwait Bay, which is heavily developed; and the south coast, which includes most of the industrial, commercial, and recreational areas. The topography of Kuwait is generally flat and decreases gradually from the southwest of the state toward the shoreline of the Arabian Gulf. The most prominent topographic area in the coastal zones of Kuwait is Jal Al-Zour. Its highest elevation reaches 145 m, and it is 60 km in length, north of Kuwait [13]. Kuwait's coastal zones occupy an area of about 1,315 km 2 , extending nearly 10 to 20 km from the shoreline. Due to the different oceanographic characterizations in coastal region, and to facilitate the interpolation processes, this study divided the coastal zones of Kuwait into five zones: northern coast (zone 1), middle coast (zone 2), southern coast (zone 3), Boubyan Island (zone 4), and Failakah Island (zone 5).

Elevation Datasets
The DEM was utilized in this research for two reasons. First, it was spatially enhanced by [11]. Second, it is highly accurate in flat areas [14], which conforms to the topography of Kuwait's coastal zone. GDEM, released in October 2011, has a spatial resolution, vertical accuracy, and horizontal accuracy of 30, 17, and 75 m, respectively. [11] matched elevation data retrieved from GDEM with 793 elevation points using regression analysis. They found that the relationship between elevations points and GDEM data was statistically significant and could be represented as the following: Elevation = -0.0012GDEM 2 + 1.2058  By applying this model to GDEM data, the nominal vertical error was greatly reduced from 17 m to 5 m. They recommended using the enhanced DEM in the absence of highly accurate data, such as LiDAR data, because it is the best available source in terms of vertical accuracy. The first elevation data used in the interpolation process was obtained from the enhanced DEM after applying this model.
The second elevation data used in the interpolation process were obtained from topographic maps. [15] produced a contour map of Kuwait with a scale of 1:400,000 but without any spatial information; thus, the vertical and horizontal accuracies were calculated based on U.S. national map accuracy standards [16]. The vertical and horizontal accuracies for this map were 10 m and 203 m, respectively. The other topographic maps used were produced by the Defense Ministry of Kuwait, with scales of 1:100,000 and 1:25,000, using in situ data and aerial photography. Both of these maps have vertical and horizontal accuracies of 0.8 m and 5 m, respectively. However, elevation points from topographic maps with scales of 1:400,000 and 1:100,000 were used to conduct interpolation.
In order to examine the vertical accuracy of the interpolated DEM, two elevation dataset collected from GPS (N= 101) and topographic maps (1:25,000; N= 48) were matched with the elevation data of the interpolated DEM in exact locations. The spatial parameters of these data are better than input data in terms of vertical and horizontal accuracies. In their research, [14] mentioned that any number between 20 and 150 reference points is more than is needed to guarantee reliable statistical accuracy measures for vertical assessment. The similarity between the interpolated DEM and measured elevations from the reference data was investigated using the linear regression coefficient (R 2 ) and was graphically illustrated.

Interpolation Processes
Both deterministic and geostatistical interpolation techniques were used to generate the continuous surface of elevation data. The deterministic interpolation techniques included IDW and RBF, whereas the kriging method was applied as a geostatistical interpolation technique.

The Inverse Distance Weighting Method
Inverse Distance Weighting (IDW) method is a wellknown interpolation technique that has been widely applied in different science fields. Inverse distance weighting predicts non-sampled locations using linearly weighted combinations of a set of sample points. This method was designed following [17] first law of geography, which states that near things are more related than distant things. Therefore, the closest points to the non-sampled location are weighted more greatly than distant points [18,19]. The method is expressed as the following equation: where z(x 0 ) is the interpolated value, n represents the total number of sample data values, x i is the ith data value, h ij is the separation distance between the interpolated value and the sample data value, and ß denotes the weighting power.

The Radial Basis Function MethodRadial Basis
Function (RBF) method covers a large series of exact, moderately quick interpolators (e.g., thin-plate spline (TPS), spline with tension (SPT), multiquadric function (MQ), completely regularized spline (CRS), and inverse multiquadratic function (IMQ)) that are used as a basic equation that is dependent on the distance between the predicted and the measured point [20]. Unlike the IDW method, the RBF predicts values above the maximum measured values and below the minimum measured values [21].

The Kriging Method
The kriging method is a sophisticated version of inverse distance interpolation that weights how much importance should be given to each neighbor location. There are several kinds of kriging, including simple kriging, ordinary kriging, and universal kriging. Ordinary kriging (OK) has been widely applied in interpolation processes due to its simplicity and prediction [22]. The weights of OK are derived from the kriging equations using a semivariance function. Semivariograms are commonly represented by plotting the difference squared between the values of each pair of locations on the y-axis relative to the distance separating each pair on the x-axis.
In this study, the following equation (2) was used to calculate the semivariogram: where z(x i ) and z(x i + h) are variable values on the x i, and x i +h points, respectively, and N(h) is the number of paired sample points separated by distance h.

Accuracy Assessment
The interpolation methods applied to generate a high spatial resolution DEM were compared using cross validation. This was conducted by removing one data location and predicting the associated data using the data of the remaining locations. The accuracy of interpolation was done by calculating the deviations of interpolated elevation values from corresponding measured values in term of the Root Mean Square Error (RMSE): where z(xi) is the observed value at location I, z*(xi) is the interpolated value at location i, and n is the sample size. Most researchers who have study DEMs and their applications proposed and used statistical measures, including the RMSE, to evaluate DEM reliability [23,24]. Moreover, the correlation coefficient (r) was applied in order to examine the relationship between the predicted elevation points from the interpolated DEM and measured values taken from input spatial datasets. Higher values of (r) indicated higher significant correlations between datasets.

Determination Errors Using Cross Validation
The spatial resolution of the interpolated DEM was defined as 5 m rather than 30 m prior to the interpolation processes. The study area was divided into five major classes: northern coast (zone 1), middle coast (zone 2), southern coast (zone 3), Boubyan Island (zone 4) and Failakah Island (zone 5). Each zone was subjected to separate interpolation analysis. In the northern coast (zone 1), the most suitable interpolation prediction with a low RMSE (1.55) was generated by the RBF (spline with tension). However, in the middle coast (zone 2), the RBF (multiquadric) achieved the best result with a low RMSE (1.56). In the southern coast (zone 3) and Boubyan Island (zone 4), the most suitable interpolated DEM was produced using the RBF (spline with tension) interpolation method, with an RMSE of 1.32 and 1.1, respectively. On the other hand, the IDW technique achieved the best result (RMSE = 1.67) compared to other interpolation methods used to generate DEMs in the Failakah Island area (zone 5).
Notably, using kriging as an interpolation technique did not demonstrate comparable results, and its estimations were mostly unrealistic. The deterministic interpolation methods reveal more accurate estimations for producing DEMs in the coastal zones of Kuwait than the geostatistical interpolation methods.
After combine the DEMs of each zone, the average error of the interpolation prediction was significantly decrease (RMSE=1.44) comparing the original spatial data used to run the interpolation. [11] reduced the vertical error of GDEM by matching remotely sensed data with ground elevation points using regression analysis. The RMSE of their regression model was 5 m. Therefore, the RMSE data in the present study add to the RMSE of the above-mentioned study to obtain the overall error of the interpolated DEM, as shown in Table 2. The overall RMSE of the interpolated DEM was 6.44, indicating an increase of 28% compared to the original data (RMSE = 5).  The degree of association between predicted and measured elevation points was examined using Pearson's correlation coefficient (Figure 3). The correlation coefficients between tested data were significant, with p values below 0.01. The elevation data of both the predicted and measured values in the northern, middle, and southern coasts were highly correlated. Correlations between predicted and measured elevation points in both the Boubyan and Failakah Islands existed but were less associated.

Visual Enhancement of DEM
The visualization of DEM was substantially improved, as shown in Figure 4. The elevation profile ( Figure 4A) in the original DEM had a spatial resolution of 30 m, which clearly shows lesser elevation bars, as in Figure 4C. In comparison, the interpolated DEM had a smaller pixel size, and the elevation bars were greater than the original DEM. The average pixel size of the original data was 30 m, but after the interpolation analysis, the pixel size was reduced to 5 m. The spatial resolution of the interpolated DEM was improved by 83% compared to the original data. The advantages of the spatial interpolation analysis amplify the effectiveness of the DEM on a local scale, in the absence of accurate elevation data.

Vertical Assessment of the Digital Elevation Model
The reliability of the interpolated DEM was subjected to vertical assessment using GCP gathered from two different sources: GPS and a topographic map (1:25,000). The elevation points of these datasets were matched with the elevation data extracted from the interpolated DEM using linear regression coefficient (R 2 ) at exact locations.
Elevation points from the interpolated DEM and GCP extracted from a topographic map ( Figure 5) and in situ elevation data ( Figure 6) were significantly correlated. In the first dataset, the elevation points from the topographic map were correlated with the interpolated DEM (R 2 = 0.88), as shown in Figure 5. However, the correlation was higher between in situ elevation data and elevation data from the interpolated DEM (R 2 = 0.94, Figure 6). These important correlations between both datasets make the interpolated DEM more reliable and effective to use in environmental fields in the coastal zone of Kuwait in absent of accurate DEM.

Conclusions and Recommendations
The interpolation process provides variety of methods that enhance the spatial resolution of DEM. the interpolation methods were carried out on different spatial datasets including original DEM, counter maps, and elevation measurements that cover the five zones of Kuwait's coastal area. The study shows that the interpolated DEM using deterministic methods such as RBF has RMSE lower than interpolated DEM using geostatistical methods. Moreover, the correlation coefficient (r) shows the significant association between predicted and measured elevation points with p values below 0.01 particularly in the zone 1, zone 2, and zone 3. The study demonstrates that elevation points from interpolated DEM matched significantly with elevation points extracted from high spatial resolution topography map and field measurements. The estimated error of interpolated DEM is about 1.44 m and that exceeds the error of spatial dataset employed in interpolation process but the spatial resolution of the interpolated DEM was enhanced significantly from 30 m to 5 m enabling further spatial analysis from national to local scale. With the lack of accurate DEM data, the interpolated DEM produced in this study may serve as a useful alternative to assess the geographic, geological and environmental issues that threaten the coastal area of Kuwait.