SRTM resample with short distance‐low nugget kriging

The shuttle radar topography mission (SRTM), was flow on the space shuttle Endeavour in February 2000, with the objective of acquiring a digital elevation model of all land between 60° north latitude and 56° south latitude, using interferometric synthetic aperture radar (InSAR) techniques. The SRTM data are distributed at horizontal resolution of 1 arc‐second (∼30 m) for areas within the USA and at 3 arc‐second (∼90 m) resolution for the rest of the world. A resolution of 90 m can be considered suitable for the small or medium‐scale analysis, but it is too coarse for more detailed purposes. One alternative is to interpolate the SRTM data at a finer resolution; it will not increase the level of detail of the original digital elevation model (DEM), but it will lead to a surface where there is the coherence of angular properties (i.e. slope, aspect) between neighbouring pixels, which is an important characteristic when dealing with terrain analysis. This work intents to show how the proper adjustment of variogram and kriging parameters, namely the nugget effect and the maximum distance within which values are used in interpolation, can be set to achieve quality results on resampling SRTM data from 3” to 1”. We present for a test area in western USA, which includes different adjustment schemes (changes in nugget effect value and in the interpolation radius) and comparisons with the original 1” model of the area, with the national elevation dataset (NED) DEMs, and with other interpolation methods (splines and inverse distance weighted (IDW)). The basic concepts for using kriging to resample terrain data are: (i) working only with the immediate neighbourhood of the predicted point, due to the high spatial correlation of the topographic surface and omnidirectional behaviour of variogram in short distances; (ii) adding a very small random variation to the coordinates of the points prior to interpolation, to avoid punctual artifacts generated by predicted points with the same location than original data points and; (iii) using a small value of nugget effect, to avoid smoothing that can obliterate terrain features. Drainages derived from the surfaces interpolated by kriging and by splines have a good agreement with streams derived from the 1” NED, with correct identification of watersheds, even though a few differences occur in the positions of some rivers in flat areas. Although the 1” surfaces resampled by kriging and splines are very similar, we consider the results produced by kriging as superior, since the spline‐interpolated surface still presented some noise and linear artifacts, which were removed by kriging.


Introduction
The Shuttle Radar Topography Mission (SRTM), a joint project of the National Aeronautics and Space Administration (NASA), the National Geospatial-Intelligence Agency (NGA-DoD) and the German Aerospace Centre (Deustches Zentrum für Luft-und Raumfahrt -DLR) was flow on Space Shuttle Endeavour in February 2000, with the objective of acquire a digital elevation 1 model of all land between 60 • north latitude and 56 • south latitude, using Interferometric Synthetic Aperture Radar (InSAR) techniques (Farr and Kobrick, 2000;van Zyl, 2001;Rabus et al., 2003;Farr et al., 2007). The data is distributed as regular-grid DEMs, at horizontal resolution of 1 arcsecond (aprox. 30 meters) for areas within the conterminous USA and at 3 arc-second (aprox. 90 meters) for the rest of the world.
SRTM data is now been widely used as source for DEMs. As a resolution of 90m can be considered suitable for small to medium-scale analysis (e.g. 1:100 000 or smaller), it is too coarse for more detailed purposes. If no data is available with higher detail, one alternative is to interpolate the SRTM data at a finer resolution. This operation will not increase the level of detail of the original DEM, but it will lead to a surface with higher coherence of angular properties (i.e., slope, aspect) between neighbouring pixels (Valeriano et al., 2006), an important characteristic when dealing with digital terrain analysis.
This work intents to show how the proper adjustment of variogram and kriging parameters (namely the nugget effect and the maximum distance within which values are used in interpolation -hence short distance-low nugget kriging) can be set to achieve quality results on resampling SRTM data from 3" to 1" resolution. Examples are presented for a test area in western USA and include different adjustment schemes (changes in nugget effect value and in the interpolation radius) as well as comparisons with the original 1", with the National Elevation Dataset DEMs, and with other interpolation methods (splines and inverse distance weighted).
The basic concepts for using kriging to resample terrain data, which will be discussed in the following sections, are: • Work only with the immediate neighbourhood of the predicted point, due the high spatial correlation of the topographic surface and omnidirectional behaviour of variogram in short distances; • Add a very small random variation to the coordinates of the points prior to interpolation, to avoid punctual artifacts generated by predicted points (pixels) with the same location than original data points; • Use a small value of nugget effect, to avoid smoothing that can obliterate terrain features.

Study area
Selection of the study area was based on the availability of 1" and 3" SRTM data -which limits the possibilities to the conterminous USA -on the morphological characteristics, such as presence of hills and flatlands, and on lack of urban features, which could introduce bias in the analysis. The test area is located in southeastern California, close to the western limit of the Mojave Block, delimited by the San Andreas and Garlock Faults (figure 1). The main geomorphological units are the southern Sierra Nevada north of the Garlock Fault, the Transverse Ranges west of the San Andreas Fault and the flatlands of the western Mojave Desert. According to Dibblee Jr. (1967), the area is composed mainly of crystalline pre-Tertiary rocks, Tertiary sedimentary and volcanic rocks and Quaternary sediments and local basaltic flows.

Dataset preparation
All kriging calculations were made using the R statistical language (Ihaka and Gentleman, 1996;Grunsky, 2002) as base system and the gstat package (Pebesma, 2004) for geostatistical functions. Since the present version of gstat cannot deal properly with data in Latitude-Longitude format, the DEM was exported as points in [x,y,z] format. These points were projected to UTM coordinate system (zone 11, Northern Hemisphere) and then imported into the R workspace for variogram modelling and kriging. GRASS-GIS (GRASS Development Team, 2007;Neteler and Mitasova, 2004) was used for raster/vector conversion and coordinate system projection.
One should avoid using a projected raster DEM as input data for resampling (instead of points in [x,y,z] format) because the coordinate projection process will create linear artifacts due pixel shifting and cell size adjustment, which can be seen in figure 2. The area depicted corresponds to the larger inset in figure 6B. The left image is a shaded relief map of the 1" SRTM DEM projected to UTM using nearest neighbours resampling, to preserve the original values, and the right image was generated by converting the Latitude-Longitude data into vector points, projecting these points to UTM, and converting them back to raster format. White areas in figure  2B correspond to "no data" pixels. It is clear that the light grey lines running in a NE-SW direction in figure 2A were introduced by the coordinate projection, in order to fill the "linear voids'. These artifacts are expected to occur systematically and to be more closely spaced with increasing latitudes.
Figure 2: Artifacts generated by the coordinate projection process. A) Shaded relief image of the 1" SRTM DEM projected to UTM. No vertical exaggeration, lighting from 315 • , 30 • above horizon. B) Elevation map created from conversion of projected vector points. White areas correspond to "no data" pixels.

Variogram modelling
The variogram is a tool that allows to describe quantitatively the variation, in space, of a regionalised phenomenon. Elevation data are usually expected to be highly spatial dependent (high similarity of data at short distances) and often have variograms with a region of low slope near the zero distance, that can be best modelled by a Gaussian model (Burrough, 1987). Variograms calculated with linear trend residues of topographic data usually have adequate fits to classical variogram models, which present a clear and defined sill (Valeriano, 2002); residues of trend surface analysis are used to guarantee geostationarity of data being modelled.
Exploratory analysis of the test area, with variograms in four directions (0 • , 45 • , 90 • and 135 • ) show that when the whole area is considered (figure 3a), the variograms have different ranges and sills, but when only the initial part of the curve is taken into account (figure 3b), the curves are very similar, which allows the use of a omnidirectional variogram for kriging, if only short distances from the origin are used.
According to Valeriano et al. (2006), the noise present in such data, that is, low similarity of data at distances close to the grid size, can be evaluated as the rate of nugget effect in variograms. This value can be calculated from the vertical error of SRTM elevation data, which has a theoretical vertical accuracy of ±16m (van Zyl, 2001;Rabus et al., 2003), and an absolute error of 9.0m for the North America (Rodriguez et al., 2006). The reported absolute error represents a 90% error in Gaussian statistics, so the standard deviation (which represents approximately 68% of the values) will be 1.69 times smaller than the 90% error (Waltham, 1995), that is, 5.32m. Since the nugget effect is a measure of the [semi]variance of data, its value will be the square of the standard deviation, that is, 28.3m 2 . This value was rounded to 30m 2 in our tests.
Due the high spacial dependency of terrain data, it is needed to limit the area within which observations are used for prediction, in order to work with small data variance and save computational cost. This can be done by limiting the maximum number of neighbouring points used by the interpolation function or by limiting the search radius of the function. In this work, we decided to constrain the number of neighbouring points, since the effect of this operation is similar to defining a square window around each point, except for the edges of the area.
Kriging honours the data values at their locations, so if a point to be estimated coincides with the location of a point in the original data, the predicted surface will have the same value as the original one at that location (Deutsch and Journel, 1992). This can create punctual artifacts in the resampled surface, where some pixels have the same value of the original DEM and the surrounding ones have interpolated values ( fig.4). Even if we consider that the 3" DEM was created by decimation of the 1" data (also known as thinning, that is, each value of the 3" DEM corresponds to the value of the central pixel in a 3x3 window of the 1" data) this situation should be avoided, because kriging is known to smooth the overall range of the data and it will also filter some of the noise present in SRTM models, so the resampled surface should be similar to the 1" surface, but not exactly equal. This problem can be solved by adding a very small random variation to the coordinates of the points prior to interpolation. One point that must be taken into account is the presence of voids in SRTM data; these voids can be successfully filled by other techniques like the Delta-Surface method of Grohman et al. (2006). One drawback of our method is that voids larger than the distance to the most distant point in the neighbourhood will not be filled, and will remain as artifacts in the resultant surface.

Results an Discussion
In order to test the influence of nugget effect and number of neighbouring points in the final results, we performed tests with nugget values of 10m 2 , 30m 2 , 50m 2 and 100m 2 , and neighbourhoods of eight (3x3 window), 24 (5x5 window), 99 (10x10 window) and 399 (20x20 window) points. The complete variogram of the area and the variogram adjusted to the initial portion of the graphic, used for interpolation, are presented in figures 5A and 5B.
The 3" data points were also gridded in GMT Smith, 1991, 1998), using continuous curvatures splines in tension (Smith and Wessel, 1990) and Inverse Distance Weighted (IDW). These methods are widely used and can be found in most GIS packages, what makes them suitable to be used for comparison with the kriging approach. To compare the interpolated surfaces with DEMs created by different methods, we selected the 1" National Elevation Dataset (NED), produced from vectorised contour lines (Gesch et al., 2002).
Examples of the analysis that were carried out are presented in figures 6 and 7. For a better visualisation of the effects of changing variogram parameters in the interpolated surface, DEM images are shown as shaded relief maps with lighting positioned at N315 • and with 30 • of inclination, no vertical exaggeration. The original 1" SRTM data is presented in figure 6A and the original 3" data in figure 6B, both with artifacts created by the coordinate projection process. The larger inset in figure 6B corresponds to the area enlarged in figures 2 and 7 and the smaller inset to the area of figure 4. The surface interpolated with kriging (30m resolution), is shown in figure 6C, and the surface interpolated with continuous curvatures splines in tension, in figure 6D. It is possible to see in the kriged surface that the coordinate projection artifacts are not present, that the noise was removed and that essentially all topographic features observed in the original 1" image can be identified. Although the surface produced by splines in tension is very similar to the kriging result, one can note the presence of noise and that the coordinate projection artifacts are still present.
In figures 7A-D, the results of changing the nugget effect values are presented. Figure 7A is a subset of the final interpolated surface with nugget=10m 2 , fig. 7B has nugget=30m 2 , fig. 7C has nugget=50m 2 and fig. 7D has nugget=100m 2 . It can be seen how the changes in this value lead to a smoother surface, with less topographical information. Figure 6: Shaded relief maps. A) original data at 1" resolution, with artifacts created by coordinate projection; B) original data at 3" resolution, larger inset is enlarged in figs.2 and 7 and smaller inset is enlarged in fig.4; C) surface interpolated with kriging at 30m resolution; D) surface interpolated with splines in tension at 30m resolution. SRTM data courtesy of NASA/NGA/USGS.

8
The effects of changing the number of neighbouring points (nmax ) is presented in figures 7E-H; fig. 7E has nmax=8, fig. 7F has nmax=24, fig. 7G  has nmax=99 and fig. 7H has nmax=399. This parameter not only influences the smoothness of the kriged surface, but also strongly determines the computational time involved, since the size of covariance matrices grows up exponentially. According to Guth (2006), SRTM 1" data compares much more closely with simulated 2" NED than with 1" NED. Following his methodology, we created a simulated 2" NED DEM by decimating the 1" data and used a plot of mean slope value versus elevation to compare the DEMs in this study. In figure 8, curves plot according to four groups, from left to right: a) original 3" SRTM (black dashed line) and 1" SRTM resampled with IDW (grey solid line); b) 1" SRTM resampled with kriging (black solid line) and with splines (grey dot-dashed line); c) simulated 2" NED (grey solid line) and original 1" SRTM (black dot-longdashed line); d) original 1" NED (black dotted line).
The original 1" SRTM data follows the observations made by Guth (2006), and plots almost atop the curve for the simulated 2" NED, while the original 1" NED plots far to right (higher mean slopes). The curve for 3" SRTM data plots far to the left (lower mean slopes), accompanied by the interpolated 1"-IDW data, which shows that inverse distance weighted resulted in an underestimation of slopes, thus it should not be used to resample DEMs to higher resolutions. Both the curves for SRTM data interpolated by kriging and by splines plot very closely, and between the curves for 3" SRTM and 1" SRTM, which shows a clear improvement of the DEM. Figure 8: Plot of mean slope versus elevation for the DEMs used in this study: 1" SRTM, 3" SRTM, 1" NED, 2" NED (simulated), 1"-kriging, 1"splines and 1"-IDW.
As a last comparison between the original and resampled DEMs, the drainage network was derived from all surfaces, using the A T least-cost search algorithm (implemented in GRASS in the r.watershed command, Ehlschlaeger, 1989). A threshold value of 100 cells was used for surfaces with 1" resolution and 33 cells for the 3" SRTM, which should represent approximately the same area. Figure 9 shows the calculated streams in two subsets of the study area, one in a region of higher slopes (same area of figures 2 and 7) and one in a region of lower slopes, where the noise present in SRTM data is more likely to induce errors in the hydrologic modelling. The "reference" drainage was derived from 1" NED (which is filtered for artifact removal) and is presented as greyscale raster backgrounds overlaid by the vectorised streams of each DEM.
Drainages derived in the hilly area (figs. 9B-E) are similar for all 1" surfaces and have good agreement with the network derived from the 1" NED, although the streams for the IDW surface are not as "clean" as the others, with duplicated and less smooth channels (central region of fig. 9E, for instance). In the low relief area, the main streams are represented in all surfaces, but not even the original 1" SRTM data was able to delineate the watersheds in the same way that the 1" NED. It is possible to see a large trunk river that was incorporated to the adjacent watershed to east (centrallower region of fig. 9G). These two watersheds were correctly identified in the surfaces resampled by kriging and spines, although some differences in the position of the main river, were some streams in the splines surface follow more closely the reference ones. Figure 9: Calculated drainage network. Grey raster backgrounds represent the "reference" streams, derived from the 1" NED, and are overlaid by vectorised versions of the streams derived for each DEM.

Conclusions
In this paper we propose the use of kriging interpolation to resample Digital Elevation Models to a higher resolution. Resample of DEMs with kriging interpolation may be a laborious task, since it involves variogram modelling prior to interpolation and care must be taken on all steps of the process, but the quality of final results were considered satisfactory.
Using of few neighbouring points allows the user to perform a good adjust of the variogram model just on the initial part of the curve, but voids larger than the distance to the most distant point in the neighbourhood will became anomaly areas or will remain unfilled; on the other hand, a larger number of points will dramatically increase computational time and produce a smoother surface, obliterating small terrain features. Nugget effect will also act as a smoothing factor; in our tests, a value of 10m 2 was sufficient to eliminate noise.
The method presented here is capable of producing surfaces with welldefined peaks and ridges, without noise, and also remove linear artifacts sometimes present in original SRTM data, but it is not suitable for void filling. In case of large voids in the data, one good approach is to first fill in the voids and then resample.
The curves of mean slope versus elevation shows that IDW resulted in an underestimation of slopes, thus it should not be used to resample DEMs to higher resolutions. The curves for SRTM data resampled by kriging and by splines plot between the curves for 3" SRTM and 1" SRTM, which shows an improvement of the DEM.
Drainages derived from the surfaces interpolated by kriging and by splines have good agreement with streams derived from the 1" NED, with correct identification of watersheds, even though a few differences occur in the positions of some rivers in flat areas.
Although the 1" surfaces resampled by kriging and splines are very similar, we consider the results produced by kriging as superior, since the splineinterpolated surface still presented some noise and linear artifacts, which were removed by kriging.