a b s t r a c tThis paper recommends computational procedures for employing auxiliary maps, such as maps of drainage patterns, land cover and remote-sensing-based indices, directly in the geostatistical modeling of topography. The methodology is based on the regressionkriging technique, as implemented in the R package gstat. The computational procedures are illustrated using a case study in the south-west part of Serbia. Two point data sets were used for geostatistical modeling: (1) 2051 elevation points were used to generate DEMs and (2) an independent error assessment data set (1020 points) was used to assess errors in the topo-DEM and the SRTM-DEM. Four auxiliary maps were used to improve generation of DEMs from point data: (1) distance to streams, (2) terrain complexity measured by standard deviation filter, (3) analytical hillshading map and (4) NDVI map derived from the Landsat image. The auxiliary predictors were significantly correlated with elevations ðadj:R 2 ¼ 0:20) and DEM errors ðadj:R 2 ¼ 0:27Þ. By including auxiliary maps in the geostatistical modeling of topography, realizations of DEMs can be generated that represent geomorphology of a terrain more accurately. In addition, downscaling of a coarse 3 arcsec SRTM DEM using auxiliary maps and regression-kriging is demonstrated using the same case study. A methodological advantage of regression-kriging, compared to splines, is the possibility to automate the data processing and incorporate multiple auxiliary predictors. The remaining open issues are computational efficiency, application of local regression-kriging algorithms and preparation of suitable auxiliary data layers for such analyses.