% Template GRASS newsletter - Article
% Language: Latex
%

% Head
\setkeys{Gin}{width=1\textwidth}
\begin{figure*}
   \centering
  \includegraphics{gnl_figures/osgeo_call}
\end{figure*}


\def\RR{\textsf{R}\xspace}


\title{Resampling SRTM 03''-data with kriging}
\label{grohmann}
%\subtitle{}
\author{by Carlos Henrique Grohmann}

\maketitle

\section{Introduction}
\label{sec:intro}
 
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 \citep{valeriano+etal: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 \footnote{\url{http://seamless.usgs.gov}} and use proper tools to import SRTM data, as pointed by \cite{neteler:2005}.
\section{Geology and geomorphology of the study area}
\label{sec:studyarea}
 

The area used as example is located in southeastern Brazil, southern region of S\~ao Paulo State (Fig.~\ref{fig:area}). 
In general terms, local geology consists of NE-SW trending metapelitic and metacalcareous rocks (Fig.~\ref{fig:geology}) of Precambrian age, deformed by the Brasiliano/Panafrican Orogenic Cycle (600-450 Ma) \citep{campanha+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 non-carbonatic rocks and lowered karstic zones, gives origin to mixed recharge systems, with important allogenic water input \citep{karmann+ferrari:2000}.

\setkeys{Gin}{width=\textwidth}
\begin{figure}[hbt]
  \centering
  %name of your graphic, without the path AND in EPS format:
  \includegraphics[width=\textwidth]{images/area}
  %caption of the figure
  \caption{Landsat 7 ETM+ image of the study area. UTM Coordinate System, Zone 22, Southern Hemisphere}
 % label of the figure, which has to correspond to \ref{}:
  \label{fig:area}
\end{figure}

\setkeys{Gin}{width=1\textwidth}
\begin{figure}[hbt]
   \centering
   %name of your graphic, without the path AND in EPS format:
   \includegraphics[width=\textwidth]{images/geology}
   %caption of the figure
   \caption{Simplified geology of study area \citep{campanha:2003}.}
   %label of the figure, which has to correspond to \ref{}:
   \label{fig:geology}
\end{figure}

\section{Variogram modelling}
\label{sec:kriging}

\textit{In this section we will deal with variogram modelling and kriging of a GRASS raster within the \RR statistical language environment. If you are not familiar with the \RR interface, please refer to \cite{bivand:2005}. Kriging procedure is partially based on GRASS-WIKI\footnote{\url{http://grass.gdf-hannover.de/wiki/How_to_interpolate_point_value_using_kriging_method_with_R_and_GRASS_6}} and on the Short draft introduction to R/GRASS interface\footnote{\url{http://www.geog.uni-hannover.de/grass/statsgrass/grass_geostats.html}}.} \\


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 \citep{valeriano+etal:2006}. Data with smooth variations (such as water levels and land-forms) often have variograms with a region of low slope near the zero distance, that can be best modelled by a Gaussian model \citep{burrough:1987}. 

According to \cite{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: {\tt karst} (the original dataset at 90m resolution), {\tt karstsub} (a subset, 90m resolution), and {\tt karst30} (same extents as {\tt karst}, 30m resolution).
\begin{smallverbatim}
>g.region n=7289760 s=7274730 w=726480 e=741870 
res=90 save=karst
>g.region n=7284000 s=7280130 w=731630 e=735860 
res=90 save=karstsub
>g.region n=7289760 s=7274730 w=726480 e=741870 
res=30 save=karst30
\end{smallverbatim}

We start working with the {\tt karstsub} region. Start the R session, call needed libraries, set up grid parameters and import the srtm raster file.
\begin{smallverbatim}
>R
>system("g.region region=karstsub");
>library(spgrass6);library(spatial);library(gstat);
>G <- gmeta6();
>grd <- GridTopology(cellcentre.offset=
c(G\$west+(G\$ewres/2), G\$south+(G\$nsres/2)), 
cellsize=c(G\$ewres, G\$nsres), 
cells.dim=c(G\$cols, G\$rows));
>mask_SG <- SpatialGridDataFrame(grd, 
data=list(k=rep(1, G\$cols*G\$rows)), 
proj4string=CRS(G\$proj4));
>srtm <- readFLOAT6sp("srtm_v2");
> names(srtm)
[1] "srtm_v2"
\end{smallverbatim}

Now we can calculate a variogram and see how it looks (Fig.~\ref{fig:varg_orig}).
\begin{smallverbatim}
>variog1<-variogram(srtm\$srtm_v2~coords[,1]+
coords[,2],loc=srtm, srtm);
>plot(variog1);
\end{smallverbatim}
%\setkeys{Gin}{width=1\textwidth}
\begin{figure}[hbt]
   \centering
   %name of your graphic, without the path AND in EPS format:
   \includegraphics[width=\textwidth]{images/variogram_original}
   %caption of the figure
   \caption{Calculated variogram}
   %label of the figure, which has to correspond to \ref{}:
   \label{fig:varg_orig}
\end{figure}

We can fit a first model "by eye"(Fig.~\ref{fig:varg_eye1}).
\begin{smallverbatim}
>vrg.eye<-(vgm(psill=11000,model="Gau",range=800,
nugget=50));
>plot(variog1, model=vrg.eye);
\end{smallverbatim}
We can check if our eyes are good and ask {\tt gstat} to calculate the model parameters.
\begin{smallverbatim}
>vrg.fit<-fit.variogram(variog1,vrg.eye);
>vrg.fit;

  model     psill    range
1   Nug  378.0942   0.0000
2   Gau 9982.4331 673.9686
\end{smallverbatim}

%\setkeys{Gin}{width=1\textwidth}
\begin{figure}[hbt]
   \centering
   %name of your graphic, without the path AND in EPS format:
   \includegraphics[width=\textwidth]{images/variogram_eye1}
   %caption of the figure
   \caption{Calculated variogram and adjusted Gaussian model}
   %label of the figure, which has to correspond to \ref{}:
   \label{fig:varg_eye1}
\end{figure}
If we want, we can adjust the parameters according to those given by {\tt 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.~\ref{fig:varg_eye2}). 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.
\begin{smallverbatim}
>vrg.eye2<-(vgm(psill=9000, model="Gau", range=550, 
nugget=5));
>plot(variog1, model=vrg.eye2);
\end{smallverbatim}
%\setkeys{Gin}{width=1\textwidth}
\begin{figure}[hbt]
   \centering
   %name of your graphic, without the path AND in EPS format:
   \includegraphics[width=\textwidth]{images/variogram_eye2a}
   %caption of the figure
   \caption{Gaussian model adjusted to the initial part of the curve}
   %label of the figure, which has to correspond to \ref{}:
   \label{fig:varg_eye2}
\end{figure}

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 {\tt gmeta6()} parameters and import the {\tt 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 {\tt jitter}, and create a new object of class {\tt SpatialPointsDataFrame}.


\begin{smallverbatim}
>system("g.region region=karst30");
>G <- gmeta6();
>grd <- GridTopology(cellcentre.offset=
c(G\$west+(G\$ewres/2), G\$south+(G\$nsres/2)), 
cellsize=c(G\$ewres, G\$nsres), 
cells.dim=c(G\$cols, G\$rows));
>mask_SG <- SpatialGridDataFrame(grd, 
data=list(k=rep(1, G\$cols*G\$rows)), 
proj4string=CRS(G\$proj4));
>srtm <- readFLOAT6sp("srtm_v2");

>coords<-coordinates(srtm);
>jcoords <- cbind(jitter(coords[,1]), 
jitter(coords[,2]));
>cat<-as.data.frame(srtm\$srtm_v2);
>srtm2<-SpatialPointsDataFrame(jcoords,cat,
proj4string=CRS(G\$proj4));
\end{smallverbatim}

Now we can (finally) use kriging to interpolate the new surface. The {\tt maxdist} parameter specifies that only points within this distance are used for interpolation. We can visualise it with {\tt image} (Fig.~\ref{fig:img_krig}).
\begin{smallverbatim}
>OK_pred <- krige(srtm\$cat~1, loc=srtm2, 
newdata=mask_SG, model=vrg.eye2, maxdist=210);
>names(OK_pred);
>image(OK_pred,"var1.pred");
>image(OK_pred,"var1.pred",xlim=c(731630,735860),
ylim=c(7280130,7284000), col=topo.colors(80))
\end{smallverbatim}
%\setkeys{Gin}{width=\textwidth}
\begin{figure}[hbt]
   \centering
   %name of your graphic, without the path AND in EPS format:
   \includegraphics[width=\textwidth]{images/image_kriging}
   %caption of the figure
   \caption{Interpolated surface, {\tt image} display}
   %label of the figure, which has to correspond to \ref{}:
   \label{fig:img_krig}
\end{figure}

The {\tt levelplot} function of the {\tt lattice} library lets you control the aspect ratio and draws a legend by default (Fig.~\ref{fig:krig_lplot}).
\begin{smallverbatim}
>library(lattice);
>levelplot(OK_pred\$var1.pred~OK_pred@coords[,1]+
OK_pred@coords[,2], OK_pred, aspect = "iso", 
main = "ordinary kriging predictions", xlab="", 
ylab="",scales = list(y = list(rot = 90)),
col.regions=topo.colors(80))
\end{smallverbatim}
Let's send the raster back to GRASS, and finish our \RR session for now. 
\begin{smallverbatim}
>writeRast6sp(OK_pred,"karst.krig",
zcol="var1.pred");
>quit()
\end{smallverbatim}
\section{Results and Discussion}
\label{sec: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 \cite{neteler:2005}.

%\setkeys{Gin}{width=\textwidth}
\begin{figure}[hbt]
   \centering
   %name of your graphic, without the path AND in EPS format:
   \includegraphics[width=\textwidth]{images/image_kriging_levelplot}
   %caption of the figure
   \caption{Interpolated surface, {\tt levelplot} display}
   %label of the figure, which has to correspond to \ref{}:
   \label{fig:krig_lplot}
\end{figure}


%\setkeys{Gin}{width=\textwidth}
\begin{figure}[hbt]
   \centering
   %name of your graphic, without the path AND in EPS format:
   \includegraphics[width=\textwidth]{images/srtm_orig}
   %caption of the figure
   \caption{Original SRTM V2 data}
   %label of the figure, which has to correspond to \ref{}:
   \label{fig:srtm_orig}
\end{figure}

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~\ref{fig:srtm_orig} is the original SRTM data, with intrinsic noise, artifacts and voids; figures \ref{fig:srtm_rst} and \ref{fig:srtm_krig} 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.

%\setkeys{Gin}{width=1\textwidth}
\begin{figure}[hbt]
   \centering
   %name of your graphic, without the path AND in EPS format:
   \includegraphics[width=\textwidth]{images/srtm_diff}
   %caption of the figure
   \caption{Difference map of krigged and RST values}
   %label of the figure, which has to correspond to \ref{}:
   \label{fig:srtm_diff}
\end{figure}

In order to compare both interpolations, the difference between them was calculated with {\tt r.mapcalc} as {\tt krigged\_map} {\it minus} {\tt RST\_map}. In Figure~\ref{fig:srtm_diff}, 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 {\tt r.neighbours}). From Figure~\ref{fig:stdev_graph}, 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 {\tt kriging\_SD} {\it minus} {\tt RST\_SD}, Figure~\ref{fig:stdev_map}),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 \ref{tbl:table1} 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 \RR environment. Also, a brief comparison with results obtained by RST was made.
%\setkeys{Gin}{width=\textwidth}

\begin{figure*}[hbt]
\begin{center}
 \begin{minipage}[t]{0.45\textwidth}  %[b] => Adjusts to \caption
  \begin{center}
   %name of your graphic, without the path AND in EPS format:
   \includegraphics[width=\textwidth]{images/srtm_rst}
   %caption of the figure
   \caption{Shaded Relief of RST interpolated DEM, illuminant at 315$^{\circ}$, 30$^{\circ}$ above the horizon}
   %label of the figure, which has to correspond to \ref{}:
   \label{fig:srtm_rst}
  \end{center}
 \end{minipage}
\hspace{.05\linewidth}                % Space between figures
 \begin{minipage}[t]{0.45\textwidth}  %[b] => Adjusts to \caption
  \begin{center}
   %name of your graphic, without the path AND in EPS format:
   \includegraphics[width=\textwidth]{images/srtm_krig}
   %caption of the figure
   \caption{Shaded Relief of kriging interpolated DEM, illuminant at 315$^{\circ}$, 30$^{\circ}$ above the horizon}
   %label of the figure, which has to correspond to \ref{}:
   \label{fig:srtm_krig}
   \end{center}
 \end{minipage}
\end{center}
\end{figure*}

\section{Conclusions}
\label{sec: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 fine-tuning 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 {\tt 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 ({\tt r.fillnulls}) and then resample with kriging.

%\setkeys{Gin}{width=1\textwidth}
\begin{figure}[hbt]
   \centering
   %name of your graphic, without the path AND in EPS format:
   \includegraphics[width=\textwidth]{images/stdev_graph}
   %caption of the figure
   \caption{Density functions for standard deviation}
   %label of the figure, which has to correspond to \ref{}:
   \label{fig:stdev_graph}
\end{figure}
%\setkeys{Gin}{width=\textwidth}
\begin{figure}[hbt]
   \centering
   %name of your graphic, without the path AND in EPS format:
   \includegraphics[width=\textwidth]{images/stdev_map}
   %caption of the figure
   \caption{Difference map for standard deviations}
   %label of the figure, which has to correspond to \ref{}:
   \label{fig:stdev_map}
\end{figure}

\begin{table*}
\begin{small}
\centering
\begin{tabular}{|c|c|c|c|c|c|c|c|c|}
\hline
%% row 1
 & krig\_DEM & krig\_slope & krig\_aspect & krig\_stdev & RST\_DEM & RST\_slope & RST\_aspect & RST\_stdev\\\hline
%% row 2
Min & 71.37 & 0 & -180 & 0.466 & 66.55 & 0 & -179.999 & 0.9144 \\\hline 
%% row 3
1st Qu. & 430.87 & 9.491 & -112.845 & 11.545 & 430.42 & 10.29 & -111.784 & 13.4169 \\\hline
%% row 4
Median & 566.56 & 15.667 & -5.742 & 17.734 & 567.06 & 16.74 & -5.602 & 19.6106 \\\hline 
%% row 5
Mean & 551.38 & 15.713 & -11.547 & 18.39 & 551.38 & 16.69 & -11.034 & 20.2738 \\\hline
%% row 6
3rd Qu. & 690.73 & 21.495 & 77.662 & 24.189 & 690.65 & 22.61 & 78.142 & 26.0236 \\\hline
%% row 7
Max & 1005.2 & 54.47 & 179.997 & 85.725 & 1005.6 & 49.58 & 179.995 & 71.6148 \\\hline
\end{tabular}
\label{tbl:table1}
\caption{Statistical summary of interpolation products.}
\end{small}
\end{table*}

\section{Acknowledgements}
\label{sec:thanks}
 

This work was supported by FAPESP grant 04/06260-5 (C.H.G. Doctor's Degree Fellowship).
The author is thankful to Roger Bivand, M\'arcio Valeriano, Markus Neteler, Hamish Bowman and Dylan Beaudette for discussions, ideas and support on the subject.


\begin{footnotesize}
\begin{thebibliography}{99}

\bibitem[\protect\citeauthoryear{Bivand}{Bivand}{2005}]{bivand:2005} 
R. Bivand (2005)
\newblock Interfacing GRASS 6 and R.  
\newblock {\em GRASS-News}, 3:11-15.

\bibitem[\protect\citeauthoryear{Burrough}{Burrough}{1987}]{burrough:1987}
P. A. Burrough (1987)
\newblock Spatial aspects of ecological data. In: Jongman, R.H.; ter Braak, C.J.F.; Tongeren. O.F.R. (Eds), 
\newblock {\em Data analysis in community and landscape ecology}. Pudoc, Wagenigen, pp.213-251.

\bibitem[\protect\citeauthoryear{Campanha and Sadowski}{Campanha and Sadowski}{1999}]{campanha+sadowski:1999}
G. A. C. Campanha and G. R. Sadowski (1999)
\newblock Tectonics of the southern portion of the Ribeira Belt (Apia\'i Domain)
\newblock {\em Precambrian Research}, 98:31-51.

\bibitem[\protect\citeauthoryear{Campanha}{Campanha}{2003}]{campanha:2003}
G. A. C. Campanha (2003)
\newblock {\em Mapa Geol\'ogico da Folha Itarar\'e (SG-22-X-B) em escala 1:250.000.}\\
\newblock \url {http://www.igc.usp.br/pessoais/ginaldo/}
\newblock \url {itarare_ld/Itarare2003.zip}

\bibitem[\protect\citeauthoryear{Karmann and Ferrari}{Karmann and Ferrari}{2000}]{karmann+ferrari:2000}
I. Karmann and J. A. Ferrari (2000)
\newblock Carste e cavernas do Parque Estadual Tur\'istico do Alto Ribeira (PETAR), sul do Estado de S\~ao Paulo. In: Schobbenhaus,C.; Campos,D.A.; Queiroz,E.T.; Winge,M.; Berbert-Born,M. (Eds.) 
\newblock {\em S\'itios Geol\'ogicos e Paleontol\'ogicos do Brasil}.\\
\newblock \url {http://www.unb.br/ig/sigep/sitio043/sitio043.htm}

\bibitem[\protect\citeauthoryear{Neteler}{Neteler}{2005}]{neteler:2005}
M. Neteler (2005)
\newblock SRTM and VMAP0 data in OGR and GRASS
\newblock {\em GRASS-News}, 3:2-6.

\bibitem[\protect\citeauthoryear{Pebesma}{Pebesma}{2004}]{pebesma:2004} 
E. J. Pebesma (2004)
\newblock Multivariable geostatistics in S: the gstat package.
\newblock {\em Computers and Geosciences}, 30: 683-691.

\bibitem[\protect\citeauthoryear{Valeriano}{Valeriano}{2002}]{valeriano:2002}
M. M. Valeriano (2002)
\newblock Modelos digitais de eleva\c c\~ao de microbacias elaborados com krigagem.
\newblock {\em Information and Documentation Service (SID), INPE}, Technical Report INPE-9364-RPQ/736, 54pp.

\bibitem[\protect\citeauthoryear{Valeriano et al.}{Valeriano et al.}{2006}]{valeriano+etal:2006}
M. M. Valeriano, T. M. Kuplich, M. Storino, B. D.  Amaral, J. N. Mendes Jr. and D. J., Lima (2006)
\newblock Modeling small watersheds in Brazilian Amazonia with shuttle radar topographic mission-90m data.
\newblock {\em Computers and Geosciences, in press, available on-line}.\\
\newblock \url {http://dx.doi.org/10.1016/j.cageo.2005.10.019}

\end{thebibliography}
\end{footnotesize}


\address{Carlos Henrique Grohmann\\
  Institute of Geosciences \\
  University of S\~ao Paulo, Brazil\\
  \url{http://www.igc.usp.br}\\
  \email{guano AT usp br}}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: main_document.tex
%%% End: 