Resampling SRTM 03”-data with kriging

SRTM is now been widely used as source for DEMs. The data is distributed at horizontal resolution of 30 meters (aprox. 1arcsec) for areas within the U.S.A. and at 90 meters resolution (aprox. 3arcsec) for the rest of the world. A resolution of 90m can be considered suitable for small or medium-scale analysis, but it is too coarse for more detailed purposes. The present alternative is to interpolate the DEM at a finer resolution. It won’t increase the level of detail of the original DEM, but it will lead to a surface where there is coherence of angular properties (i.e., slope, aspect) between neighbouring pixels (Valeriano et al., 2006), an important characteristic when dealing with terrain analysis. The purpose of this article is to present the steps necessary to improve the resolution of a DEM using variogram modelling and kriging, as well as a brief comparison of the results with those obtained with interpolation by Regularised Splines with Tension (RST). GRASS users must be aware of the half-pixel shift in SRTM "finished" data from USGS web site 4 and use proper tools to import SRTM data, as pointed by Neteler (2005).


Introduction
SRTM is now been widely used as source for DEMs. The data is distributed at horizontal resolution of 30 meters (aprox. 1arcsec) for areas within the U.S.A. and at 90 meters resolution (aprox. 3arcsec) for the rest of the world.
A resolution of 90m can be considered suitable for small or medium-scale analysis, but it is too coarse for more detailed purposes. The present alternative is to interpolate the DEM at a finer resolution. It won't increase the level of detail of the original DEM, but it will lead to a surface where there is coherence of angular properties (i.e., slope, aspect) between neighbouring pixels (Valeriano et al., 2006), an important characteristic when dealing with terrain analysis.
The purpose of this article is to present the steps necessary to improve the resolution of a DEM using variogram modelling and kriging, as well as a brief comparison of the results with those obtained with interpolation by Regularised Splines with Tension (RST).
GRASS users must be aware of the half-pixel shift in SRTM "finished" data from USGS web site 4 and use proper tools to import SRTM data, as pointed by Neteler (2005).

Geology and geomorphology of the study area
The area used as example is located in southeastern Brazil, southern region of São Paulo State (Fig. 1). In general terms, local geology consists of NE-SW trending metapelitic and metacalcareous rocks (Fig. 2) of Precambrian age, deformed by the Brasiliano/Panafrican Orogenic Cycle (600-450 Ma) (Campanha and Sadowski, 1999) and affected by Cretaceous brittle tectonics and basic dike emplacement.
Karstic landscapes developed over the carbonatic rocks, with altimetric differences up to 700m between non-carbonatic (pelitic, psamitic and granitic) crests and karstic valley bottoms; the structural pattern of the area, alternating elongated ranges of noncarbonatic rocks and lowered karstic zones, gives origin to mixed recharge systems, with important allogenic water input (Karmann and Ferrari, 2000).

Variogram modelling
In this section we will deal with variogram modelling and kriging of a GRASS raster within the R statistical language environment. If you are not familiar with the R interface, please refer to Bivand (2005). Kriging procedure is partially based on GRASS-WIKI 5 and on the Short draft introduction to R/GRASS interface 6 .
The variogram is a tool that allows to describe quantitatively the variation, in space, of a regional phenomenon. Elevation data are usually expected to be high spatially dependent (high similarity of data at short distances); 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 (Valeriano et al., 2006). Data with smooth variations (such as water levels and landforms) often have variograms with a region of low slope near the zero distance, that can be best modelled by a Gaussian model (Burrough, 1987).
According to Valeriano (2002), variograms calculated with linear trend residues of topographic data have adequate fits to classical variogram models, which present a clear and defined sill. Residues of trend surface analysis are used to guarantee geostationarity of data being modelled.
Since DEMs datasets can be very large, variogram calculation and kriging interpolation can become a very time-consuming tasks. The basic idea is to choose a representative subset of the study area, and calculate the variogram over this subset. This variogram model will then be used by kriging to interpolate the whole dataset at a finer resolution.
In this example, I worked with three regions: karst (the original dataset at 90m resolution), karstsub (a subset, 90m resolution), and karst30 (same extents as karst, 30m resolution).

Figure 3: Calculated variogram
We can fit a first model "by eye" (Fig. 4).
>vrg.eye<-(vgm(psill=11000,model="Gau",range=800, nugget=50)); >plot(variog1, model=vrg.eye); We can check if our eyes are good and ask gstat to calculate the model parameters.  If we want, we can adjust the parameters according to those given by gstat. Since the initial part of the curve of the Gaussian model and the nugget effect will both have the effect of smoothing the interpolated surface, we can concentrate on adjust just this part of the model (Fig. 5). Using a nugget effect value of zero would result in a very noisy surface, while a large value produces a smooth surface, maybe at the cost of loosing detail.
>vrg.eye2<-(vgm(psill=9000, model="Gau", range=550, nugget=5)); >plot(variog1, model=vrg.eye2); Figure 5: Gaussian model adjusted to the initial part of the curve After the variogram model is adjusted, it is time to interpolate the new DEM. First, we need to set the region back to the full area of the DEM, but with finer resolution, re-set (in this case, overwrite) the gmeta6() parameters and import the srtm_v2 raster. Then, we need to add a very small random variation into the coordinates of the points, to avoid interpolation artifacts later. We add the noise with jitter, and create a new object of class SpatialPointsDataFrame.

Results and Discussion
Interpolation by RST has been successfully used for void filling of SRTM data, and can also be used to interpolate the DEM at a finer resolution. For more on this subject, please refer to Neteler (2005). As a means of comparison with the results obtained with kriging, the SRTM dataset was also interpolated at 30m resolution using RST (with default values). Figure 8 is the original SRTM data, with intrinsic noise, artifacts and voids; figures 10 and 11 are shaded relief maps for the interpolated DEMs with RST and kriging, respectively.
From the figures above, we can see that RST map has a sharper look and does a better job about void filling (note the dark areas in the upper right corner of kriging map), but linear artifacts are still present. The krigged map has a smoother appearance, with crests not so well defined as in the RST map, but the linear artifacts were essentially removed. In order to compare both interpolations, the difference between them was calculated with r.mapcalc as krigged_map minus RST_map. In Figure 9, positive values (areas where kriging values are higher than RST ones) have colors ranging from blue to red and negative values have colors from cyan to yellow; The greatest differences are in valley bottoms and in crests due the smoothing behaviour of the kriging function.
To better evaluate the differences of interpolated products, standard deviation (SD) maps where calculated with a 7x7 moving window (using r.neighbours). From Figure 12, we see that the SD of kriging interpolation is a little smaller than the SD of RST. If a difference map is made from the SD map (as kriging_SD minus RST_SD, Figure 13),the presence of artifacts in RST interpolation is enhanced; the extreme values are related to areas of voids in the original data, that were not completely filled by kriging. Table 10 presents a statistical summary of the interpolated DEMs and their derivatives slope, aspect and standard deviation (all derivatives were calculated with a 7x7 window).
The purpose of this article was to present a simple guide on how to increase SRTM data resolution using kriging within the R environment. Also, a brief comparison with results obtained by RST was made.

Conclusions
RST interpolation is a strong method, suitable for void filling, and produces good results, although it does not eliminates the artifacts inherent to SRTM data. A more detailed study concerning the finetuning of RST parameters may point interesting directions on the use of this tool. Kriging interpolation may be a more laborious task, since it involves variogram modelling prior to interpolation. Care must be taken on all steps of the process. The use of the maxdist option allows the user to perform a good adjust of the variogram model just on the initial part of the curve, but larger voids will became anomaly areas or will remain unfilled. Nugget effect will act as a smoothing factor; a small value is sufficient to eliminate noise, while a larger one may obliterate terrain features.
In case of large voids in SRTM data, one possible approach is to first fill the voids with RST (r.fillnulls) and then resample with kriging.