Geografisk Tidsskrift, Bind 99 (1999)
Generation of a digital elevation model of Mols Bjerge, Denmark, and georeferencing of airborne scanner data
Anne Jacobsen, Niels Drewes, Michael Stjernholm & Thomas Balstrøm
The study compares various digital elevation models (DEM) and georeferencing methods for the purpose of orthorectifying airborne scanner data. The study was performed using topographical maps on the scale of 1:25,000 produced by the Danish National Survey and Cadastre and a Compact Airborne Spectrographic Imager (casi) image. The study area. was a topographically complex area in Denmark, Mols Bjerge. The generation ofthe digital elevation model was performed using a discretised thin plate spline technique modifled to allow thefitted DEM tofollow abrupt changes in the terrain. The accuracy ofthe DEM was primarily evaluated by comparing the original and the modelled contour lines. The georeferencing was performed using tegiunul mudifitd vtrsiuns of polynumial dislortion modelling andtriangle warping, applying bothground controlpoints and DEM features for panoramic and topographical correction in the process. It was concluded that the maps had an accuracy thai allow - ed for a DEM with a grid resolution of lOm and that triangle warping performed betler than the polynomial distortion modelling in the case of airborne scanner data since the image distortions åre related to local distortions superimposed on medium to regional scale distortions. Triangular irregular network(TlN) resampling performed betler than regional modelling especiallly in hilly areas outside the nadir line where extracted tie points described the landscape features. The overall accuracy ofthe georeferencing was 3-4 pixels (6-8m), which was well within the reference map accuracy of 12-15m.
grid resolutions, topographic correction, polynomial
Jacobsen: National Environmental Research Institute,
Balstrøm: Institute of Geography, University of
Stjernholm: National Environmental Research Institute,
Danish Journal ofGeography 99: 35-45,1999
Acquisitions of airborne remote sensing data covering Denmark have increased within the last 5-6 years and there is a rising demand for accurate georeferencing of the image data. Airborne remote sensing data is generally recorded at low altitudes, using a sensor with a wide field of view (FOV), which imposes pronounced geometric distortions on the image data. The geometric distortions vary with topographical differences, even though the Danish landscape has relatively low elevations ranging from-5m to +173m. The final result of remote sensing analysis often necessitates a combination with other georeferenced data, e.g. in a geographical information system (GIS) and thus,
of geometric image distortions becomes
Remote sensing data have for a long time been georeferenced using polynomial transformation models based on corresponding pairs of ground control points (GCPs) found in the image and reference maps. Since satellite orbits åre fairly stable and image lines åre acquired at a speed of approximately 300 Hz, a transformation model in the form of a Ist or 2nd order polynomium usually works well (Schowengerdt, 1997). In cases where e.g. SPOT HRV images have been acquired in the slant mode, over rugged terrain, correction methods have been applied to
remove the effect of the terrain distortion. One method is to manipulate the GCPs in order to project the image coordinates onto a selected reference level, calculate the transformation model and the resampling grid, and then project the resampling grid onto a digital elevation model of the area.
Since 1985 a digital elevation model (DEM) of Denmark with a spatial resolution of 50m, basedon maps on the scale of 1:50,000, has been available (Frederiksen, 1985). The use of this DEM with satellite data over a study area on Fyn, Denmark, proved that image georeferencing and co-registration were improved using 3-dimensional GCPs (Sandholt et al., 1995). However, the 50m grid DEM is too coarse for some environmental applications (Larsen et al., 1999) and also for the georeferencing of airborne scanner data. For this purpose a DEM with better spatial and altimetric resolution is needed. The 2.5m contours and extreme points from the 1:25,000 maps åre respected as being the bestbasic source when creating elevation models for Denmark (Larsen et al., 1999) and these maps were used as the basis for the DEM in the study. The data quality of the maps and the accuracy of DEMs with different grid sizes åre discussed in the paper.
The acquisition of airborne low altitude data is performed with a rate of 10-30 lines per second under unstable conditions due to attitude (roll, yaw and pitch) variation caused by aircraft motion, wind speed, and turbulence. This results in image rotation and also in along and across track displacement, and therefore each image line should, ideally, be treated individually with information on its absolute orientation. The Compact Airborne Spectrographic Imager (casi) image was acquired with a system not fully optimized for the acquisition of high precision attitude data. Thus the georeferencing section of the paper has two parts: 1) evaluation of georeferencing accuracy, using polynomial resampling models correcting for elevation and panoramic effects across track basedon DEM of different grid sizes, and 2) a modification of piecewise triangle warping using GCPs and TIN vertices extracted from the grid based DEM as tie points.
The study area comprises Mols Bjerge, which is a test area within the Danish Multisensor Airborne Campaign (DANMAC) (Groomet al., 1997). Mols Bjerge is located on the Djursland peninsula in the eastern part of Jutland (figure 1). Mols Bjerge is one of the most hilly areas in Denmark with elevations up to 137 m. Mols Bjerge was formed during the last ice age 16,000 years ago. The area was located between two glacier tongues and is characterized by terminal moraines pushed from both sides forming rather long, parallel ridges in the landscape (Larsen & Kronborg 1994). Behind the terminal moraines a belt of dead-ice was situated; kames, kettles and deep erosion valleys were formed by the melt water (Landskabsplan for Molsområdet 1978). Today, a few agricultural fields, open grasslands with some scrubs and thickets, deciduous forests and coniferous plantations dominate the area (figure 1).
Due to the glacial history that can be so easily recognized in the landscape and because of the rich flora and fauna in the area, the northern part of Mols Bjerge was declared a protected area in 1977. The southern part was added to the protected area in 1992.
The DEM is developed from the basic topographical map series on the scale of 1:25,000 produced by the Danish National Survey and Cadastre (Kort og Matrikelstyrelsen, KMS). These maps include 2.5m equidistant contour lines,
locations of local depressions and peaks and information on land cover and infrastructure. The contours were originally created from photogrammetry in 1984. The general horizontal precision of the 2.5m contours is approximately 12-15m (0.5-o.6mm on the map sheets) and the vertical accuracy is reported to be close to 0.75-1 m. The location and precision of local peaks and depressions åre assumed to fall within this range. As the contour density for Mols Bjerge is very high, it is likely that the quality of the vertical accuracy is greater herethan for the rest of Denmark and maybe close to 0.5m.
Image and reference data
In the summer of 1997, a Compact Airborne Spectrographic Imager (casi) was flown over Mols Bjerge. Six images were acquired close to solår noon on the lOth of June. The casi scanner used had a 42° FOV and an Ilband spectral configuration from 400 to 900 nm. The spatial resolution was 2 m by 2 m (nadir line pixel size). The image used in the study covered approximately l km by 3 km of the central part of Mols Bjerge (figure 1). The scan line was influenced by the aircraft attitude (roll, pitch and yaw) but as the casi scanner was mounted with a gyro, roll correction was possible. The sensor was not mounted on a three axes stabilized platform so correction for aircraft pitch and yaw was not possible. GCPs for georeferencing the image data were found using the Danish Digital Orthophotos (DDO) (Kampsax/Geoplan, 1995). Differences in vegetation structure, mtersections between grasslands and forests, gravel roads and nature paths were recognizable features in the image and DDO. Control and evaluation of the georeferencing was basedon differential GPS readings using Trimble AgGPSI32 with an accuracy of +/- 2 m.
Digital elevation modelling
Several interpolation techniques have been developed and Integrated into geographical information systems during the last decades in order to produce a DEM from digitized contours. Depending upon the nature of the data source one may, for example, choose to select a minimum curvature or a kriging algorithm to produce a cell based DEM. Occasionally, a TIN-based interpolation approach may be
preferable (Burrough and McDonnell, 1998). The data from the study area showed a lack of a high number of elevation indications for local peaks and depressions. A TIN-based interpolation was, therefore, not an option as too many local depressions and peaks would have been cut off at the lowest/highest lying contour levels, causing several of the landscape peaks and valleys to be 'flattened out'.
In order to generate a DEM, based mainly on digitized contours, an interpolation procedure named TOPOGRID, mcluded m the ARC/INFO ver. 7.2.1 software from ESRI (ESRI, 1998), proved to give the most satisfying DEM for Mols Bjerge when carefully tuned. The TOPOGRID module(Hutchinson, 1988;Hutchinson, 1996)is basedon an interpolation procedure designed to take advantage of known characteristics of the landscape. TOPOGRID uses an iterative finite difference interpolation technique optimized in order to have the computational efficiency of 'local' interpolation methods such as inverse distance weighted interpolation, without losing the surface continuity of global interpolation methods such as kriging and splines. It is essentially a discretised thin plate spline technique, where the roughness penalty has been modified to allow the fitted DEM to fol lo w abrupt changes in the terrain such as streams and lakes. Therefore, it becomes possible to include the following data sources into the interpolation: contours, peaks, depressions and locations of streams and rivers (the latter two entered with twodimensional extents, only). Three parameters can then be adjusted to control the output (ESRI 1998):
1) How much drainage should be hindered given a specific threshold value (e.g. half the contour equidistance). Data points which block drainage by no more than this value åre removed.
2) The amount of error inherent in the process of converting point, line, and polygon elevation data into a regularly spaced grid - i.e. a horizontal component. The program (depending on the local slope at each data point and the grid cell size) scales this value. Large values will cause data smoothing resulting in a more generalized output grid. Smaller values will cause less data smoothing resulting in a sharper output grid more likely to contain spurious depressions and peaks.
3) The amount of
non-systematic (random) error in the
z-values of the
See below for a
general description of the tuning of the
The regional approach is a modification of the commonly used polynomial distortion models. The polynomial distortion models åre approximations ofthe coordinate transformation between reference map and input image coordinates estimated from a set of GCPs. The polynomial order is often between one and three. A linear polynomium (Ist order) works well if the input image has been processed for systematic distortions, 2nd order models work well for most remote sensing satellite problems, where the terrain relief is low and the FOV small (Schowengerdt, 1997). The casi has a large FOV and a high terrain relief relative to the flying altitude and 2nd and 3rd order distortion models åre evaluated keeping in mind that there åre a number of disadvantages using the polynomial method: 1) The polynomial transformation can only deal with a limited number of smooth changes in image geometry; 2) the method does not guarantee a correct transformation at GCP locations; 3) GCPs from locally distorted parts ofthe image may lead to inaccuracies over a wider area (Devereux, 1990).
The resampling procedure was supplemented by a resampling grid, which specified the sampling coordinate in the input image at a regular interval of the output image. During resampling, the sampling coordinates between grid points were interpolated.
The regional approach applied panoramic and topographical corrections to the ground control points and the resampling grid. The panoramic and topographical corrections were applied for each GCP pair (Ix, Iy, Rx, Ry, Rz) by projecting the input sample Ix to a reference elevation (Rref) using:
where Ix, Iy denotes the input image coordinates (sample and line), Rx, Ry, Rz denotes the reference map coordinates (sample, line and elevation), Øx is the local view angle, IFOV is the instantaneous field of view, Dx is displacement across track, pe is the local pixel width and Ix proj is the projected sample coordinate at the reference elevation (figure 2).
After calculation ofthe polynomial resampling model and the resampling grid, each grid sampling point was projected from the reference elevation to the local topographic surface using elevations read from the DEMs. The projection was done using the formulas la and le and modifying Ib to:
resampling was performed using a standard
The local approach is a modification of piecewise triangle warping (TIN resampling) where GCPs åre used to define a net of triangles. During the resampling, the calculation of the sampling position is performed by bilinear interpolation using the vertices of a given triangle as the tie point. The triangular network can be generated by e.g. Delaunay triangulation (Schowengerdt, 1997). Devereaux et al. (1990) used a modification where the initial GCPs åre supplemented by iteratively interpolated new tie points found by triangulation in both map and image space. The iteration stopped when each triangle was below a dermed area threshold. This modification does not model the panoramic and topographic distortion.
We applied the
correction for panoramic and topographic
1) For each tie
point (GCP point), the input sample coordinate
2) The projected
tie points were used to divide the
defined by the tie points into a net of
triangle were modeled as supplementary tie
triangle warping procedure.
Digital elevation modelling
The DEM is basedon two digital maps on the scale of 1:25,000; 1314 I NV Helgenæs and 1315 II SV Rønde. The conversion process from digital maps to digital elevation model is illustrated in the flow diagram figure 3.
Contour lines and 4 map corners were acquired as l bit tiff files and the contour lines were traced using ARCSCAN software (ESRI, 1998). The 5m contours were represented by solid lines and therefore easily traced. The 2.5m contours were dashed lines and particular attention was paid to ensure proper snapping of data (on the slopes). The 4 map corners were manually digitized. The contour lines were then auto-labeled using a macro program for the ARC/INFO software called AUTOKONTUR (Balstrøm and Haffner, 1997). Afterwards the two resulting coverages were georeferenced to the UTM coordinate system using 8 depressions and the 4 map corners as GCPs on each coverage. Comparison between the georeferenced traced contour lines and the original contour lines showed that the two data files fitted each map sheet well, with only few meters displacement (within the accuracy of the used maps in the x,y-plane) in some areas in the eastern part of 1315 II SV.
The two coverages were joined using the two lower map corners of 1315 II SV Rønde and the two upper map corners of 13141 NV Helgenæs as tic links. The final contour coverage covering both map sheets was then compared with the digital maps. The contour lines on 1314 I NV Helgenæs perfectly overlaid the original data, but the contours were displaced about lOm to the east and 6m to the south on the 131511 SV Rønde map. This displacement is due to the disagreement between the position of the map corners, which åre not identical on the digital maps. However, this offset is within the accuracy of the topographic maps.
elevation points, lakes, streams, the shoreline
Georeferencing of casi data using the regional approach was basedon the flow diagram illustrated in figure 4. The GCPs were sampled with image x,y-coordinates, UTM coordinates and elevation was read from the DEMs in different resolutions. The GCPs were projected onto a reference level (Om) and the geometric distortion model and grid were calculated. Then the grid points were then projected
reference level onto the DEM and the
Application of the local approach using steps 1-4 described in the theory chapter implies that vertices of a TIN characterize the terrain. DEMs of the area in the TIN form were created for each of the grid resolutions mentioned above. The approach was basedon iterative Delaunay triangulation of the vertices given by the previous iteration followed by an interpolation of the resulting TIN. Each triangle was then tested to see whether the difference between the regular grid DEM and the TIN interpolated DEM was above an acceptance level of lm. If the difference exceeded the acceptable level, a new vertex was added to the positions where the maximum divergence in both the positive and negative directions was detected. The TIN resampling of the image was done using the ENVI 3.1 software (RSI, 1998).
Results and discussions
Output from TOPOGRID
Several DEMs with different spatial resolutions and settings for various options were produced using TOPOGRID. All interpolated DEMs were evaluated by a visual comparison between generated 2.5m contours from the individual cell based DEMs and original input contours. If the input and output data overlapped, it was assumed that the interpolation had been successful. As mentioned in the methods section the input data consisted of digitized contours, peaks and depressions (in total 378 points, very few compared to the approximately 490,000 input points from the contour vertices), lakes, streams and a shore line, which acted as part of the boundary. TOPOGRID's stream, lake and enforcement options were chosen early in the experiments to promote the locations of the landscape's lowest lying areas, because the DEMs were better modelled using the locations of these features as an additional gridding feature.
As the input contours åre located very close to each other in several areas (less than lOm apart) the generated output grids generated the best contours when the spatial resolution was refined (figure 5). In general, the location of local peaks or depressions was better modelled the finer the grid size.
227 depressions were added as input to the interpolation. However, as shown in table l, the interpolation produced a much higher number of depressions even though the spatial resolution was set to a low value. Whether these
depressions åre real or false depressions within the landscape, is very difficult tojudge without directly comparing the results to the actual on site topography but in this type of landscape it is, of course, impossible to put all occurring depressions onto maps on the scale 1:25,000. So, it is possible that the number of generated sinks from the TOPOGRID runs åre close to the real conditions. Meanwhile it was found that lowering the tolerance parameter that controls the block drainage, decreased the number of generated depressions.
The number of depressions was adjusted by setting the horizontal tolerance and number of iterations (table 1). The recommended horizontal tolerance (ESRI, 1998) is half the interval of the contour lines, which in this case is 1.25 m, but as the contours åre very densely located it was decided to vary this parameter within the range 0.125-1.25 m. An increase of the horizontal error smoothened the surface and decreased the number of depressions. A decrease of the horizontal error raised the number of depressions and produced a more rough grid. Generally, a higher value than 2 did not decrease the number of depressions (a value of l is default). The vertical standard error was set to the default value of O, in all models, assuming that the input data were not expected to have any significant random vertical errors.
As an example it can be mentioned that using lOm grid resolution and a tolerance of 1.25m the enforcement routine removed approximately 17% of the 1553 generated sinks. A tolerance setting of 0.75m gave the best result leaving 995 false sinks in the model - a reduction of an additional 23%. A further reduction of 18% was achieved by changing the horizontal standard error from l to 2 thereby increasing the number of iterations to 50 and reducing the number of sinks to 740. The number of itera
tions did not have a large effect on the number of sinks beyond the values listed in table 1. A higher number of iterations resulted in only 1-2% fewer false sinks located outside the boundary polygon. Increasing the number of iterations smoothed the surface but only in areas where the data density was low. An example of a DEM is shown in figure 6, see next page.
A better DEM for Mols Bjerge would be appreciated, but according to Larsen et al. (1999) a second generation DEM for the whole of Denmark will not be available before the year 2001. This new model will be basedon revised 2.5m contours and additional features such as streams and lakes with attached elevation values. This product will form a much better basis for the production of a more precise DEM for Mols Bjerge, than presented in this article. However, it is doubtful whether the Danish National Survey and Cadastre will recommend the producing of a DEM with a greater resolution than lOm basedon this forthcoming data set.
The evaluation of the accuracy of the regional georeferencing and the effect of different DEM gridsizes is based on the goodness of fit (residual mean square, RMS, table 2) for the polynomial models. Table 2 shows the RMS for different DEM gridsizes and different polynomium orders using 45 GCPs.
Table 2 shows that the regional approach using 2nd order polynomium modelling and raw GCPs, the traditional approach for satellite image georeferencing, has been improved by l pixel using 3rd order polynomial modelling and GCPs modified for displacement across track due to
scan angle and topography. The RMS is within the same pixel accuracy regardless of DEM grid resolution and distortion model. The RMS is 2-4 pixels (6-8m), which is well within the accuracy of the maps and orthophotos. The overall correlation between the results shown in table 2 is most likely due to the faet that the information provided to the DEM generation has an accuracy in the order of 10-15 meters.
The RMS only
provides a part of the information on the
quality of the
georeferencing process since the georefer
encing, after the
calculation and polynomial modelling,
includes the calculation and manipulation of the resampling grid. In this context it is important that the resampling grid oversamples the variation in topography. To ensure this, a grid spacing of 10 pixels (30m) was used. The resampled image, using a 3rd order polynomium model, 3m
DEM and nearest
neighbour resampling is shown in figure
The local approach is dependent on GCPs to describe the attitude of the aircraft/sensor and a high number of tie points to describe the terrain. We therefore combined 65 GCPs with TIN vertices (6500 points) extracted from the regular DEM. TIN resampling basedon a 3m DEM is shown in figure 7, subsets l and 3. Figure 7 shows 4 subsets; Subsets l and 2 åre a relatively flat area along the nadir line, including road crossings used as GCPs. Subsets 3 and 4 show a very hilly area. No GCPs were placed in this area. The crosses represent the gravel roads and foot paths in the area tracked with DGPS and superimposed on the image. The dropouts åre due to inadequate readings of radio beacons for differential positiomng. The f irs t two subsets in figure 7 show that the regional and local approach have largely the same results along the nadir line and in areas where GCPs have been sampled. Subsets 3 and 4 show that the TIN resampling improves georeferencing considerably further away from the nadir line and in hilly areas where no GCPs have been used for the modelling. Since there were no GCPs found in the area, the closeness in correpondence between the tracked paths and the TIN resampling is due to the inclusion of vertices extracted from the DEM.
Polynomial models cannot, in general, be expected to georeference airborne scanner data satisfactorily. Airborne scanner data åre subject to local distortions superimposed on medium to regional frequency distortions. TIN resampling, using automated selection of tie points from a DEM, identified local landscape features and therefore performed better than 3rd polynomial modelling. Further TIN resampling has the benefit of improving poor correspondence
in one area by inserting a new point in the area. Adding another pair of GCPs will also improve the resampled result using the regional approach. This will, however, most likely move the area of weak correspondence te another region. The drawback of TIN resampling is that only part of the image, inside the convex hull of tie points, can be resampled.
The results of the digital elevation modelling showed that the best DEM outputs were generated when the output grid resolution was set to 10 or even sm. The production of a 3m DEM proved useful for georeferencing purposes but a DEM with a smaller grid size than lOm is generally not recommended since the original contours from the 1: 25,000 maps have a horizontal precision of ± 12-15m. As the input contours åre very densely located in Mols Bjerge, it is assumable that the positional horizontal error is in the vicinity of ± l Om. The vertical accuracy of the input data was reported to be close to 0.75-I.OOm. As it is possible to regenerate the 2.5m contours quite precisely from the best generated DEMs, it is assumed that the average vertical error of the generated DEM is still within this range.
Georeferencing of the casi scanner data using digital elevation data did not add any significant improvement in accuracy using the regional approach. The precision was 2-4 pixels independent of DEM grid size. The accuracy of the distortion models was due to the overall precision of the maps and orthophotos. The improvement of RMS by l pixel using a higher order polynomium and georeferenced GCPs is not regarded to be significant. Including automatically extracted vertices from a DEM and local resampling using a straightforward Delaunay triangulation turned out to be substantially important for the overall georeferencing accuracy. Comparison between the resampled images using 3rd order polynomium modelling and triangulation proved that georeferencing is especially improved in off nadir hilly areas without GCPs.
It is concluded that is possible to produce a DEM for a very undulating landscape with a grid spacing of lOm with reasonable horizontal and vertical confidences using maps on the scale 1:25,000 and that TIN georeferencing is successful and superior to polynomial distortion modelling when tie points have been automatically extracted from a DEM at even higher resolution.
The authors would like to thank The Danish Space Board for financing the DANMAC program, Assoc. Prof. Birger Hansen for support and encouragement during the digital elevation modelling process and Ib Krag Petersen for field assistance and preprocessing of the DGPS data.
Balstrøm, T. & Haffner, J. (1997). AUTOKONTUR - an ArcTool-AML for the preparation of contours as an input source to digital elevation models. Proc. 12* ESRI European UserConference, Sept. 29-Oct. l, 1997. Copenhagen. Publ. on CD-ROM, Informi GIS, Denmark.
& McDonnell, R.A. (1998): Principles of
Fuller, R.M., Carter, L., & Parsell, R.J. (1990):
(1985): Generating DEMs from contours using
Groom, G., Mikkelsen, H., Skriver, H., Hansen, 8., Thompson, A., & Nielsen, A. A. (1997): Activities and Early Results from DANMAC: The DANish Multisensor Airborne Campaign. Proc. Third International Airborne Remote Sensing Conference and Exhibition, 7-10 July 1997, Copenhagen, Denmark.
Hutchinson, M.F. (1988): Calculation of hydrologically sound digital elevation models. Third International Symposium on Spatial Data Handling, Sydney. Columbus, Ohio: International Geographical Union.
Hutchinson, M. F.
(1996): A locally adaptive approach to the
Jacobsen, A., Carstensen A.R., & Kamper, J. (1993): Mappingof Satellite Derived Surface Albedo on the Midtlaugkat Glacier, Eastern Greenland Using a Digital Elevation Model and SPOT HRV Data. Geografisk Tidsskrift 93: 6-18.
Jacobsen, A., Heidebrecht, K.B. & Nielsen, A.A. (1998): Monitoring Grasslands Using Convex Geometry and Partial Unmixing -a Case Study. In: Proc. lsl EARSeI Workshop on Imaging Spectroscopy. Eds. Shaepman, M. et al., 6-8 October 1998, Ziirich, Switzerland.
(1995): Danmarks Digitale Ortofoto (DDO).
Molsområdet (1978): Landskabsplan for Molsområdet.
Balstrøm, T. & Jacobi, O. (1999): Towards a second
Larsen, G. &
Kronborg, C. (1994): Det mellemste Jylland. En
områder af national geologisk interesse.
RSI (1998): ENVI
3.1. Research Systems, Inc., and Better Solutions
Sandholt, 1., Fog, 8., Poulsen, J.N., Stjernholm, M., & Skriver, H. (1995): Classification of agricultural crops in Denmark using ERS-1 SAR and SPOT imagery. In: Sensors and Environmental Applications of Remote Sensing, Proc. 14th EARSeL Symp., Goteborg, Sweden. Balkema, Rotterdam, pp. 37-44.
R.A. (1997): Remote Sensing, models and