REMOVAL OF TREE OFFSETS FROM SRTM AND OTHER DIGITAL SURFACE MODELS

The recently completed 1 second Digital Elevation M odels (DEMs) for Australia are based on the 1 secon d Shuttle Radar Topographic Mission (SRTM) elevation data. The SRTM d ata was corrected by removing voids, striping, tree offsets and random noise and finally by integrating mapped drainage li n s. This paper describes the removal of the tree o ffs ts, which was a crucial step in the production of a credible bare-earth elevatio n model and was one of the most technically challen ging aspects of the project, and the possible application of the methods to other di gital surface model (DSM) sources. Methods for the removal of tree offsets rely on map s of tree presence/absence from sources such as rem otely-sensed imagery, and the height offsets are computed from the DEM at the bou ndaries of tree patches. Tree offsets over most of Australia were successfully removed, but were underestimated in areas of extens ive forest cover and poorly estimated where the map ping of tree patches did not match the patterns of offsets in the SRTM elevations . The tree offset removal methods could be applied to the near-global SRTM DSM to produce a near-global b are-earth product, provided that a suitable map of tree presence or de nsity can be compiled from satellite remote sensing a d other sources. The process could be improved by using supplementary tree-heigh t information from ICESat or other sources. High resolution global DEMs other than SRTM are beco ming available, notably ASTER GDEM and TANDEM-X. Both t ose products are subject to offsets due to vegetation i n the same way as SRTM. The tree offset removal meth ods developed for SRTM could be adapted to the characteristics of these an d other DSMs to provide a largely automated process ing ystem to derive bareearth DEMs from new sources.


INTRODUCTION
The release of the Shuttle Radar Topographic Mission (SRTM; Farr et al., 2007) data for Australia in 2005 heralded a significant step forward in resolution, detail and consistency for elevation models over that continent.Previous continental DEMs had been prepared from topographic mapping at scales of 1:2.5M to 1:100k with the most recent being the Geodata 9 second DEM version 3.0 (http://ga.gov.au/topographicmapping/digital-elevation-data.html).The impact of this improvement in resolution is enormous.Most hillslopes are on the order of 200 m long so changing from 250 m to 30 m resolution brings the ability to resolve the basic ridge-valley structure of the land surface that is crucial to understanding the flow of water and materials through the landscape and everything that follows from that, in particular the patterns of soil depth and properties and the distribution of vegetation communities.
While the improved resolution provided obvious benefits, the defects in the SRTM data were also readily apparent: in contrast to the smooth, complete and hydrologically enforced 9 second DEM the SRTM data was noisy, had missing data (voids), contained stripes, was subject to significant offsets in vegetated areas and did not capture river channels.In short, it is a noisy digital surface model (DSM) rather than a bare-earth model (DEM or DTM).This led to a collaborative project between CSIRO, the Bureau of Meteorology, Australian National University and Geoscience Australia to deal with each of those defects extending from 2005 through 2011.Access to the 1 second SRTM data (which is not routinely released for territories outside the USA) was provided by DIGO, Australia's defence mapping agency.
Widespread diagonal striping with a wavelength of about 800 m was detected and treated using a 2-dimensional Fourier transform in a custom-built tool that allowed manual identification of the striping frequency and orientation in ¼ × ¼ degree tiles (Read et al., in prep).Voids were removed using a modification of the Delta Surface Fill method (Grohman et al., 2006), with the 9 second Geodata DEM providing infill data (Read et al., in prep).Tree offsets were removed using the methods described in this paper.Random noise was reduced using an adaptive smoothing method (Gallant, 2011) that smooths to a greater or lesser degree in response to the noise amplitude (estimated from the SRTM DSM) and the local relief.Hydrological connectivity was enforced with a modified version of the ANUDEM program (Hutchinson, 2011) using mapped 1:250k stream lines modified to fit the SRTM DSM in higher relief areas (Dowling et al., in prep).The resulting set of four products (a cleaned DSM, bare-earth DEM, smoothed DEM-S and hydrologically enforced DEM-H) are now publicly available through Geoscience Australia's National Elevation Data Framework Portal (http://nedf.ga.gov.au) at both 1 second and 3 second resolutions (except for the DSM, which is only available for government use at the 1 second resolution, and DEM-H, which is only available at 1 second resolution due to the difficulty of generalising while retaining the necessary hydrological connectivity).This paper describes the methods developed to remove tree offsets from the SRTM DSM, which was a crucial step in the production of a credible bare-earth DEM for Australia.These tree offsets are clearly visible in the SRTM DSM in agricultural lands and managed forests where there are abrupt transitions between cleared and tree-covered areas.Riparian forests and woodlands are also problematic: many of Australia's most significant rivers run through very low relief terrain and are typically bordered by tall trees while the surrounding landscape is either cleared for agriculture or is naturally devoid of dense tree cover.The rivers therefore appear in the SRTM DSM as raised ridges rather than channels.Collectively these various tree offsets prevent effective use of the DSM for most analytical purposes.
Manual, or manually assisted, workflows for removing offsets from DSMs exist within commercial organisations but were considered too expensive to use across an entire continent.Fully automated methods were therefore developed to remove the offsets due to trees.The methods produced effective results over most of the vegetated parts of the continent (about half of the area of Australia required some removal of tree offsets).Some manual effort was required to select the best tree cover map for use in different areas and to evaluate the results but the processing itself is fully automated.Offsets due to buildings and other structures were not treated.
The opportunity now exists for these methods to be applied elsewhere, potentially producing a global bare-earth DEM suitable for most analytical applications.The estimated tree height offset is also likely to be valuable as a high-resolution surrogate for biomass.With suitable adaptation, the methods could also be adapted to other DSMs such as TANDEM-X to fully automate the production of bare-earth DEMs.

METHODS
The existence of tree offsets in the SRTM data was an expected characteristic of the data, and several studies have investigated the relationship between tree height and the induced offset in SRTM elevations where a reference DEM was available (e.g.Kellndorfer et al., 2004;Walker et al., 2007).At the time this work was commenced there were no published examples of estimation and removal of the offset where a reference bareearth DEM was not available, and even now the authors are aware of only one other similar effort (Köthe and Bock, 2009).
The method described here relies on an independent map of tree cover so that the areas requiring treatment can be defined.This map must be compatible with the DSM resolution and represent the conditions at the time of the DSM acquisition (February 2000 in the case of SRTM).Land cover maps derived from Landsat TM at about 30 m resolution are the obvious choice, and in Australia we were fortunate to have several such maps derived from Landsat.These maps were produced using different processes and from imagery at different dates so a manual selection process was applied to choose the best representation of the visible SRTM offsets in different parts of the country.
The main steps in the tree offset removal process are: 1. Adjustment of the tree cover patches to match edges in the SRTM DSM; 2. Characterisation of the SRTM DSM to patch edges; 3. Estimation of tree offset at patch edges; 4. Interpolation of offset to the interior of patches; and 5. Subtraction of the tree offset surface from the DSM.
A full description of the processing steps and the compilation of the tree cover map will be provided in Gallant et al (in prep).

Adjustment of the tree cover map
The correct estimation of tree offset relies on accurately identifying the location of the transition from non-tree-covered to tree-covered terrain.Mis-registration between the DSM and the tree cover map will result in under-estimates of the tree offset and ineffective removal of the offsets, so a process was developed to adjust the edges of patches in the tree cover map to optimise the match between the map and the DSM.
Several modified versions of the tree cover map were created, and the effectiveness of each version was assessed using an F ratio in a neighbourhood around each grid cell that measured the difference in average height between the tree-covered and non-tree-covered area, relative to the variation in height within both areas.For each cell, the version with the highest F ratio was chosen as the best representation of the patch edge.If the best F ratio was smaller than a minimum value, this was deemed to indicate that there was insufficient information to reliably adjust the patch edges and the original map was used; this mostly occurred in higher relief areas where the terrain variation masks the offset due to vegetation.The range of adjustment was limited to two grid cells to prevent unreliable adjustments, so the method cannot correct for large errors in the tree cover map.

Characterisation of edge response
The response of the SRTM DSM to changes in tree cover is not sharp but transitions smoothly over a distance of 3-4 grid cells (around 100 m).This smooth transition must be accounted for in the correction of offsets to avoid artefacts around tree cover patches.
The response to edges was assessed over a large (>10 6 ) sample of clearly defined edges in flat terrain and modelled using a Gaussian smoothing kernel with a length scale of 1.4 grid cells.The adjusted tree cover map was smoothed using this kernel for use in the remainder of the processing.

Estimation of tree offset at patch edges
The tree offset was estimated near every patch edge by least squares estimation based on a model of local height variations: • the estimated variance of h ˆ is less than 3 m 2 ; • h ˆ is more than twice its estimated standard deviation, since low values of h ˆ are usually caused by a poorly matching tree map; and • h ˆ is less than 25 m to exclude unreasonably large estimates that typically occur where tree patterns coincide with terrain features.The estimated parameters 3 0 â a K are not used in the subsequent process but are included in the model so that the estimated tree offset is not overly influenced by correlations between the orientation of the patch edge and the natural gradient of the land.

Interpolation of tree offset across patches
The height offset estimates obtained from the least squares estimation are only available for patch edges whereas heights are required throughout the patches, so the height estimates must be interpolated through the patches.This was achieved using an early version of the multiscale adaptive smoothing method (Gallant, 2011) that uses the variance information to weight the estimates.This produces a continuous surface of tree height offset both within and outside patches i.e. for both treecovered and non-tree-covered areas.

Subtraction of offset from DSM to produce DEM
The interpolated tree offset is multiplied by the smoothed tree cover map to produce an estimate of the tree offset that can be subtracted from the DSM to produce the bare-earth DEM:

Computational details
Most of the processing steps were implemented as AML macros within ESRI Arc Workstation.The least squares fitting was implemented as a C++ program.
Due to the size of the SRTM DSM for Australia (about 40 GB of single-precision floating point data) the processing was performed on 1×1 degree tiles using sufficient overlaps in each step that the results were consistent across tile boundaries.This tiled processing approach also allowed for distributed processing.At various stages of the project we used an 80-node processing cluster, a multi-core server and a Condor distributed processing system using idle desktop PCs.
The tree offset estimation was performed on the DSM product after removal of stripes but without voids filled, so that the interpolated data in voids would not influence height estimates.After the subtraction of estimated tree offsets, voids were filled and water bodies re-flattened as described in Read et al., (in prep).

RESULTS AND DISCUSSION
Figure 1 shows an example of the removal of tree offsets in an area where the method worked very effectively.The DSM (centre panel of Figure 1) shows marked elevation offsets that would seriously interfere with land surface attributes such as slope and with flow paths computed from the DSM.After treatment (DEM, right panel of Figure 1) there are very few artefacts due to vegetation patterns.Both the regular patterns of plantation forestry and the more erratic patterns of natural vegetation in the valley (running from the eastern to southern border of the area) are well represented in the tree cover map and the height offset is accurately modelled by the algorithms.
Figure 2 shows the estimated tree height offset for the area shown in Figure 1.The variation in estimated offset around the patch edges where the offsets are measured contrasts with the smooth variation within the patches where offsets are interpolated from the edges.The lower estimated height within the patches is due to variations in heights around the edges and is probably an under-estimate: it is more likely that trees are at least as high within the core of the patch compared to the edges.The impact of these errors is relatively minor and the method adequately deals with the most troublesome aspect of the artefacts, the abrupt changes in height at patch edges.
The techniques described here were applied to 813 1°×1° tiles covering Australia.Effective removal of vegetation offsets, assessed by visual examination of DEM, was achieved for about 90% of the vegetated area of the Australian continent.The remaining areas contain offsets that were untreated or only partially removed, or areas where the offsets were overestimated.Most of these defects are related to the tree cover mapping and have several causes including: • the mapping could not distinguish between trees and other ground cover; in some cases areas of trees are missed while in other cases areas not covered by trees are classified as tree-covered • there were significant changes of tree cover in the time between the imagery capture for the mapping and the SRTM radar mission Some of the defects are due to the algorithms, particularly overestimation of tree height offset where vegetation and terrain patterns coincide (such as tree covered rises in cleared agricultural land) or in areas of low vegetation height where the constraints on acceptable height estimates may have produced biased estimates by omitting small but correct height offsets.
Coincidence of vegetation and terrain patterns can also result in under-estimated offsets, for example in riparian forests filling inset floodplains.The interpolation of height offsets across large areas of continuously wooded landscape also appears to have resulted in under-estimation of tree height offset in much of the interior of the patches.
Figure 3 shows an area where some areas affected by tree offsets have not been effectively removed because of differences between the tree cover map and the patterns of offsets in the SRTM DSM, presumably due to changes in tree cover between the Landsat and SRTM acquisition dates.The regular patchy nature of the tree cover is due to plantation forest harvesting and replanting.Figure 3 also illustrates over-estimation of tree height offsets around the river.The river is incised, bordered by cliffs nearly 20 m high, and the algorithm has been unable to separate the co-incident effects of sudden terrain change and tree cover change thus interpreting the cliffs as tree heights and substantially over-estimating the tree height offset.The removal of this apparent offset has largely removed the cliffs and produced land heights beside the river comparable to the river surface height.
Improvements in the quality of the derived bare-earth DEM could readily be achieved by correcting parts of the tree cover map that did not correspond well with the offsets in the SRTM DSM.It is also likely that the algorithms could be enhanced to provide more accurate results.

Applications within Australia
The bare-earth DEM produced from the SRTM DSM is still quite noisy and lacks the hydrological connectivity needed by some applications.The DEM for Australia has been further processed to produce a smoothed (DEM-S) product and a hydrologically enforced product (DEM-H).From these two products a suite of commonly used terrain attributes have been derived, which are now being used for ecological, hydrological and geomorphological purposes.The stream networks and catchments derived from DEM-H will underpin the ongoing refinement of the Australian Hydrological Geospatial Fabric  (AHGF, http://www.bom.gov.au/water/geofabric/index.shtml, commonly known as the 'Geofabric'), an authoritative and rich representation of Australia's hydrological features.
The tree height offset map could be used as the basis of a high resolution map of forest height and biomass, provided suitably robust relationships can be found between height offset and actual tree height.Note that the tree height offset derived from the SRTM data is not a measure of true tree height, primarily because the SRTM radar signal scatters from the woody structure within the canopy rather than the upper margin of the canopy (Kellndorfer et al., 2004).Preliminary study of the estimated offsets suggests the estimated heights are a smaller fraction of true height in trees with a conical shape (typically conifers) than in trees with a broad canopy (such as most Australian eucalypts), presumably due to greater penetration of the radar signal into the former canopy type.

Prospects for processing of the global SRTM DSM
The processing method for estimating and removing tree height offsets in the SRTM DSM could in principle be applied to the entire SRTM dataset.The most significant obstacle to attempting this is the need for a consistent global land cover map compatible with the SRTM data in both resolution and acquisition date.Global cover maps currently exist at resolutions of about 300 m (GlobCover) and 500 m (MODIS Land Cover).Finer resolution products from Landsat might also be possible, and it is also possible that the raw SRTM radar product would contain useful information for detecting vegetation cover boundaries.
Independent measures of tree height (suitably modified to account for SRTM penetration into the canopy) could be incorporated into the height offset estimation process.This would be particularly helpful in extensively forested areas where there are likely to be spatial variations in tree height that cannot be estimated from patch edges.Either direct measurements from instruments such as the GLAS laser altimeter aboard ICESat or indirect measurements based on multispectral imagery could be used.Some promising progress in combining those two approaches was reported by Lefsky (2010) using GLAS and MODIS data to produce a global tree height map.While this is at a relatively coarse resolution it may provide sufficient information to support improved offset estimates in extensive forest areas; further development of the methods could also provide improved information.Combining the various data sources under a model-data fusion approach may also be possible, yielding a tree map, tree heights and bareearth DEM from a single process.

Application to other DSM types
While some of the aspects of the algorithm are specific to SRTM, notably the matching of response to patch edges, much of the method is in principle directly applicable to any DSM that responds predictably to vegetation cover.A bare-earth DEM is usually produced from a DSM by manually driven editing and hence tends to be quite expensive.A fully automated method, even if it is not as accurate as a manually driven method, is attractive where very large areas of DSM need to be processed, as with global or continental DSMs like SRTM.
Two obvious candidates for processing using the method presented here are the ASTER GDEM and TANDEM-X products (although the ASTER GDEM response to trees does not seem as consistent as the SRTM so may be less amenable to treatment).Both are essentially DSMs and would benefit from removal of offsets due to tree cover.

CONCLUSIONS
Removing artefacts from Digital Surface Models is a prerequisite for deriving high quality bare-earth Digital Elevation Models.The methods described here provide that capability, drawing on an appropriately fine resolution vegetation (tree) map as supporting data.The resulting bareearth DEM for Australia is now being used for a variety of ecological, hydrological and geomorphological applications.
The methods described here could be applied to the entire SRTM near-global DSM, provided a suitable tree cover map can be compiled.This would be a valuable step in realising the full value of this high resolution, high quality product and would complement other efforts such as those by CGIAR to provide a void-filled version of the SRTM DSM.
The method may also be suitable for removing offsets due to tree cover in other remotely sensed DSMs such as ASTER GDEM and TANDEM-X.Some adaptation to the different characteristics of these DSMs would be required.
h are fitted to values of DSM z and G m from a circular window of 5 cell radius around each target cell using linear least squares fitting.The estimated tree height h ˆ is retained if: • the model fits well enough (defined by 200 2 < χ );

Figure 2 .
Figure 2.Estimated tree height offset for the area shown in Figure 1.