III. Developing a Framework for Comparison
This section serves two goals. First, it indicates the methodology
by which this study compares the elevations of non-coincident
points. An overview of the bilinear interpolation method is provided,
and the implementation of this algorithm for the datasets used
in this study is indicated. The precision of this implementation
is difficult to measure, and this topic is granted some room.
The chapter discusses the statistic by which the two data sets
are initially compared. Secondly, the section shows how one set
of elevation points may be transformed horizontally or vertically
prior to interpolation. This linear transformation can be used
to attempt to improve the fit between two coincident elevation
sets and to understand the nature of systematic error in one or
another of the DEMs.
Figure 3.1. Distribution of DEM spot elevations across a
portion of the Santa Barbara DEM. For the purposes of display
within a GIS, the elevation data sets and difference maps were
converted to raster format. This involved resampling of the point
data in several cases. No analysis was done using the resultant
resampled rasters; they are employed solely for data visualization.
Once data are in a comparable format, how can they be compared?
From Figure 3.1, it is clear that the two sets of points are not
coincident in space. The 3 arc second locations are in a gridded
pattern with a slightly different orientation - and clearly different
spacing - than the 30m data. I decided to evaluate the two data
sets' characterization of the terrain at the 3 arc second point
locations. Some interpolation method must be used in order to
establish an elevation value for each of those locations from
the 30 meter data set. The use of gridded DEM structure implies
nothing about the value of locations between the known points,
unlike a TIN (Kumler, 1992), so the interpolation technique is
not defined by the data model. I chose bilinear interpolation
for two reasons. First of all, it is a straightforward method;
Li (1993) notes that "linear surfaces are usually the least
misleading" when comparing them with higher order surfaces.
It is simpler to implement and faster to compute than other methods
(e.g. thin plate splines). Since the interpolation would have
to be computed over 22,000 times for each data set, and the intention
was to run it multiple times on a relatively large number of data
sets, speed was an issue. The general format of the data in Figure
3.1, with each 3 arc second value lying within a nearly square
quadrilateral defined by the surrounding 4 30m locations, suggests
that this method will be computationally appropriate.
The general bilinear method is graphically illustrated in Figure
3.2; terms in the following discussion are identified there. The
four corner points Z1-Z4 are known elevation values, while elevation
at the point (u,v) is what we wish to calculate. First, elevation
is linearly interpolated for Z' and
Z'' at locations (u,v')
and (u,v''), respectively. The following
equations accomplish this:
(3.1)
(3.2)
The division term in both equations is the fraction of the distance
along the x1-x2 line segment that u occurs. This is used
to identify the height, based on a linear relationship in elevation
between the two known points. A final linear interpolation occurs
along the line segment u. In this case, the elevation at
location v along this line is of interest, and can be calculated
using the following equation:
(3.3)
The implementation of the formula uses some simplifying forms
to increase the speed of the algorithm, but the result is the
same. This method portrays the surface within the four points
as a curved form; Kumler (1992) likened the surface to a twisted
ladder, with sides and rungs remaining straight but not resting
in the same plane.
Figure 3.2. Bilinear interpolation. Points Z1-Z4 have known
elevations. Z(u,v) is the elevation to be estimated.
The core of the implementation was, for each of the roughly 22,000
3 arc second locations falling within a given USGS 7.5' quadrangle,
to find the 4 bounding 30m points and to perform a bilinear interpolation
for that location. The difference between this interpolated value
for the 30m characterization of terrain could then be compared
to the original 3 arc second elevation for every 3 arc second
location in the quadrangle. I was unable to find a pre-developed
algorithm that would accomplish this. The C code I wrote to implement
the comparison deserves some discussion. Since the algorithm had
to identify the bounding points for thousands of 3 arc second
locations, a data structure was required for the 30m data set
that could implement fast and efficient searches on large spatial
data sets. The structure that I settled upon, a double nested
linked list, was conceptually very similar to the DEM storage
structure used by the USGS (compare Figure 3.3 with Figure 2.3).
Figure 3.3. Linked list data structure used for the grid
search. Zx indicate both an individual elevation sample and its
corresponding easting and northing.
The double linked list models the transects of the original 30m
storage structure. Each element in the list includes an easting
value representative of values for the entire transect. Following
the datum transformation to WGS 72, the 30m data set did not maintain
exactly the same easting within each transect; this value varied
by as much as .03 meters. Each of the transect elements also includes
three links; one is to the transect element to the east, one back
to the west, and one to the first point element in the transect.
This structure is called a double linked list because all elements
have pointers to elements on both sides, rather than just to the
following element. The elevation point values themselves are stored
in nested linked lists branching off of each of the transect elements.
Each of these location elements included the easting, northing,
and elevation for that location, as well as a pointer to the next
element in the list. This structure is able to locate values very
quickly. The code maintains a pointer to the transect that was
last visited by the preceding search. The next search -- for
the bounding four 30m points of next 3 arc second location - will
begin in this transect. This speeds the search since each location
is quite probably in the same (or possibly an adjoining) transect
as the previous location.
Once the program finds the bounding four 30m points for a 3 arc
second x-y location, the elevations and coordinates of these points
are used to estimate an elevation value using the bilinear interpolation
technique discussed earlier. The difference between the 3 arc
second value and this interpolated value is calculated, and the
next 3 arc second location is input for processing. The C program
can accomplish this quickly. As an example, constructing the linked
list structure for the 176,294 points in the Cuyama Peak 30m DEM
takes about 4 seconds on Whizbang, a Dec Alpha in the department.
Processing the 27,268 3 arc second locations from the Los Angeles-w
DEM that lie near or within this quadrangle - finding the bounding
box for each point, calculating the bilinear estimate, and generating
a root mean difference fit also takes about 4 seconds on this
computer. On Yak, an IBM RISC/6000, times were 38 and 19 seconds,
respectively. On Whizbang processing was fast enough to allow
a large number of DEMs to be processed repeatedly. Fifteen DEMs
were processed to determine an estimate for the number of links
traversed by the code. For these files, this number ranged from
10,732,000 to 16,263,000, with 9 of the 15 falling above 16,200,000
links traversed.
This research used root mean square error (RMSE) to quickly evaluate
the quality of the relationship between the two DEMs. RMSE is
a common statistic for quantifying elevation error (Shearer, 1990),
and is used by the USGS as a data quality standard for many of
its products, including DEMs (USGS, 1993b.). The following formula
is a standard presentation of the RMSE:
(3.4)
where v is the is the error at each point, and n
is the number of points being evaluated. For this research, the
"error" component of the RMSE is actually a difference
value. The comparison program I wrote calculated RMSE as:
(3.5)
where a is the 3 arc second value and b is the bilinearly
interpolated value for each location. Since the number of points
evaluated is quite large (a typical quadrangle contains about
22,000 3 arc second locations), the difference in the denominator
has little impact on the RMSE statistic.
While much of the RMSE is related to differences in the two DEMs,
some fraction of it may be the result of processing error. This
error accrues from the accumulation of small rounding errors,
glitches in the source elevation data, and the ever-present possibility
of programming bugs. It is important to develop an idea of the
magnitude of this fraction. One way to do this is through the
use of small synthetic data sets. These sets are samples of planar
surfaces for which the elevation at any point can be easily calculated.
By perturbing the points in one dataset, I could generate RMSE
statistics that could also be calculated by hand. This result
is matched with the RMSE calculated by the program; any difference
indicates a problem with some aspect of the code. The program
(the current version, at least!) passes these tests with ease,
indicating that the basic logic is correct.
A second test is to compare a real-world DEM against itself.
By using a large number of actual DEM elevations and coordinates,
this test assesses the code in a more realistic operating environment.
Because the values and locations of all points are the same, the
RMSE should be 0. However, this was not the case for any of the
30m DEMs that I tested. The Santa Barbara level 1 DEM, for example,
had an RMSE of 0.1224 meters across its 176,676 bounded points.
Differences ranged from -1.735 meters to 2.024 meters, with 90%
of the differences falling between -0.007 and 0.008 m.. The Ashton
level 2 DEM, with 171,632 bounded values, had an RMSE of 0.0217.
Actual differences ranged from -1.08 m. to 1.11 m., but 90% of
the differences fell between -0.01 and 0.011 m.. None of the 30m
DEMs tested against themselves had an RMSE greater than .128 meters.
These figures are presented to provide a sense of the "slop"
in the calculation of RMSE for this research. In a sense, they
indicate the limits to the precision of this method's calculation
of difference.
There is debate about the utility of the RMSE in geographic data
assessment. This statistic measures only absolute difference at
each point, which raises two distinct questions for this application.
First of all, absolute elevation error may be immaterial as a
measure of accuracy. For example, if the entire elevation data
set is uniformly offset from the actual surface, it has retained
the characteristics of the elevation surface and terrain attributes
such as slope, aspect, and relative height are unaffected. More
broadly, this point-based measure of difference ignores any spatial
patterns exhibited by the difference surface. Difference residuals
turn out to be spatially autocorrelated for the elevation datasets
I examined. Due to the simplicity of calculation and its availability
in data specifications, however, many studies employ it as a measure
of accuracy for terrain studies (Hunter & Goodchild, 1995;
Adkins & Merry, 1994; Polidori & Chorowicz, 1993; Fisher,
1991; case studies in Shearer, 1990), and in the evaluation of
spatial data generally (Garbrecht & Starks, 1995; USGS, 1993b).
This work uses RMSE as a first cut at the problem of assessing
the relationship between the two collocated data sets. Spatial
measures were used to evaluate the distribution of the resulting
difference surfaces.
The code allows the user to transform the coordinates of the 3
arc second locations prior to searching for the bounding 30m points.
This allows the user to explore the possibility that horizontal
shifts in the data are responsible for some portion of the difference
between the two files. By interactively manipulating the transformation
and assessing the impact of the shift upon the RMSE value, the
fit can be improved. The following paragraphs will provide some
detail on this portion of the research.
The linear function developed by this code uses 6 coefficients to transform the coordinates for each point (u,v) to locations (u',v'):
3.6
3.7
Employing the coefficients allows one to shift the dataset by
a constant distance in u and/or v, and also provides
for a linear shift of each point by a distance relative to its
location in u and/or v. For example, to maintain
the original location, coefficients b and f are
set to 1 while the others are set to 0. To shift the values 30
meters to the north, d is set to 30. To make the horizontal
shift a function of the point's east-west position in the region,
b is set to some fraction. By manipulating these 6 values,
one can investigate the possibility that horizontal error is contributing
to the differences between the two datasets. A seventh coefficient,
g, is employed to change the elevation, if one wishes to
assess the possibility of constant vertical error in the comparison.
The software developed for this research queries the user for
these 7 parameters prior to doing any processing. Points are transformed,
and then the bounding four points are identified. The order is
important, since transformation may well cause 3 arc second points
to shift to locations outside of their initial bounding points,
as Figure 3.4 demonstrates.
Figure 3.4. A sample transformation which results in a different
set of bounding 30m points. Circles represent 30m locations; the
X marks the initial, pre-transformation 3 arc second location;
the boxed X indicates the location following transformation.
This section has presented important technical aspects for the
implementation of this research. Discussion of tool development
is neglected in most papers I've read, but I felt it was important
to include here. This feeling is partially due to the effort and
time I have invested in tool production, but also because I think
the decisions made in the development stage profoundly affect
the results of any experiment. Technical methodology is more than
simply background information. The tools which implement this
framework will be exercised in the following sections, and I believe
understanding the framework for comparison will aid in assessing
the findings of this research.