50K and 5K scale topomaps and aerial photo have been obtained from the Croatian State Geodetic Department (http://www.dgu.hr). Orthorectified photo was produced following the methodology explained in Rossiter & Hengl (2002). From the orthophoto we digitized on-screen land cover polygon map using the following classes: agricultural fields, fish ponds, natural forest, pasture and grassland, and urban areas. From the stereo-pairs we interpreted the generic landforms and then created a polygon map of geoforms (see also Rossiter & Hengl (2002)). Nine landform elements were recognised: summit, hill shoulder, escarpment, colluvium, hillslope, valley bottom, glacis (sloping), high terrace (tread) and low terrace (tread). From topomaps, we extracted contours and streams and water bodies. In the case of 1:50K the equidistance was 20 m in hilland and 5 m in plain, and for the 1:5K the equidistance was 5 m in hilland and 1 m in plain. From the 1:5K contours and geodetic points, the 5 m DEM has been derived using the ANUDEM (topogrid) procedure in ArcInfo and then resampled to the 25 m gird. The 30 m SRTM DEM (15'x15' block) was ordered from the German Aerospace Agency (http://eoweb.dlr.de), then resampled to the 25 grid so it can be compared with the DEM25m. IMPORTANT NOTE: According to a licence agreement, the SRTM dataset can not be distributed or used for commercial purposes outside this project.

The 5020 field measurements of elevation (land survey) collected in the 1960s by the 'Meetkundige Dienst Rijkswaterstaat'. This was used to generate the 25 m DEM25TOPO. The 5 m LiDAR-based DEM distributed by the Ministry of transportation and water management (measurements in cm). This dataset is also known as "Actueel Hoogtebestand Nederland" (AHN ) (van Heerd et al. 2008).

All input raster maps are in ArcInfo *.asc format, and the point data (tables) are in a *.dbf format. All coordinates are in the official German coordinate system, zone 3 (germany3): Transverse Mercator Projection, central meridian is 9°, false easting 3500000, Bessel 1841 ellipsoid with Potsdam datum. The bounding coordinates of the study area are: XMIN=3570000, YMIN=5708000, XMAX=3580000, YMAX=5718000. The input raster maps are available in two grid resolutions: 25 m (fine) and 100 m (coarse). The sand, silt and clay values have been determined by using the so-called texture by hand method: a surveyor distinguishes to which of the 32 texture classes a soil samples belongs to, and then estimates the content of fractions. E.g. texture class St2 has 10% clay, 25% silt and 65% sand.

The data set was obtained from the USGS National Map seamless server (http://seamless.usgs.gov). The map of soil mapping units was obtained from the Natural Resources Conservation Service (NRCS) Soil Data Mart (http://soildatamart.nrcs.usda.gov). The scripts used to predict soil mapping units and extract landforms are available via the authors website. The elevations range from 1400 to 1800 meters. There are six soil mapping units: (1) Holland family, 35 to 65% slopes; (2) Chaix-chawanakee family-rock outcrop complex; (3) Chaix family, deep, 5 to 25% slopes; (4) Chaix family, deep, 15 to 45% slopes, (5) Holland family, 5 to 65% slopes, valleys; (6) Chaix-chawanakee families-rock outcrop complex, hilltops.

All other data sets are from lidar, but the coverage varies, some data sets cover only beach and foredune area and do not include the main dune. Date of the surveys is in the name of of the file in the form YYYY_MM_DD.

Digitized from a topographic map by Ross Ihaka. These data should not be regarded as accurate.

> data(volcano)
> library(spatstat)
> LLC <- data.frame(E=174.761345, N=-36.879784)
> coordinates(LLC) <- ~E+N
> proj4string(LLC) <- CRS("+proj=longlat +datum=WGS84")
> LLC.NZGD49 <- spTransform(LLC, CRS("+init=epsg:27200"))
> volcano.r <- as.im(list(x=seq(from=2667405, length.out=61, by=10),
+     y=seq(from=6478705, length.out=87, by=10), 
+     z=t(volcano)[61:1,]))
> volcano.sp <- as(volcano.r, "SpatialGridDataFrame")
> proj4string(volcano.sp) <- CRS("+init=epsg:27200")
# str(volcano.sp)
# spplot(volcano.sp, at=seq(min(volcano.sp$v), max(volcano.sp$v),5),
+    col.regions=topo.colors(45))
> write.asciigrid(volcano.sp, "volcano_maungawhau.asc", na.value=-1)

The original topo-map DEM was produced by digitizing contour layers from two adjacent sheets of the 1:5000 topographic maps with contour interval of 5 m. Two sheets were scanned by ANATech Evolution scanner with 400 DPI resolution, then georeferenced to the Gauss-Krüger coordinate system (7th zone) and converted to a point map using a semi-automated digitalization of contour lines (I/GEOVEC, Intergraph program module). The faults obtained during automated digitalization were removed by 3D editing of contour lines. The final DEM was produced using the ArcGIS 3D analyst: first a TIN was produced, which was then converted to a regular grid of 30~m resolution. This will be referred to as the topo-DEM in further text. The original points extracted from the topo-map (51,847 points) were sub-sampled (for computational efficiency) to 2051 points.

A set of 1020 photogrammetric control points was provided with the help of the Geodetic governmental authority of Serbia. These were obtained throughout the orthophoto map production of the scale 1:1000 for local municipality. Aerial images were obtained using the RMK 21/23 analog camera with calibrated focal length f=207.96 mm. The average flying altitude was 1040 m. A stereoscopic model was produced using the WILD A10 analog stereo-restitution instrument and MapSoft2000 program package. The land-surface points were measured manually from the stereoscopic model with average lag of 15 m. This gave a total number of 46,021 points that were sub-sampled to 1020 points for faster processing. The estimated height accuracy of control points is 15 cm which allows us to use it as ground truth for the topo-DEM.