III. Kriging and Error Analysis
The Kriging interpolation technique was used to determine DEM error. As
stated previously this is usually done with ground control points, but
these are not available for this site. I would highly suggest finding a
geostatistics expert to help with the kriging process. Kriging is a
complicated interpolation procedure and an expert in this field would
be ideal. For this research, the kriging analysis was
processed in Surfer 7.0, geoeas 2.1, and ArcGIS 9.0 Geostatistical
Analysis. Point kriging can solely be done in ArcGIS, but block kriging
needs to be determined in Surfer 7.0. The image below describes
each step of the kriging process. I will then go into the detail and
background of each step.
Overall Kriging/Error Analysis Steps:
Kriging Steps:
1. First the DEM
difference is determined outside of the volcanic region. A mask
is used to clip the area outside of the volcanic region. To do this,
first create a new shapefile which will be a polygon. The steps for
creating a new shapefile and editing it are found here. Add the newly created shapefile into ArcMap. First I added a field in the attribute table (shown below).
Ultimately I wanted one shapefile to show both the mask and the outside
region. Here I used the value field to distinguish the two and will
apply that to a raster. Once the field is added, the editor
toolbar is used to edit the new shapefile. For this study I just used a
box to
identify the outline of the DEM, which is available on the Advanced
Editing Toolbar. Once the outline of
the DEM is finished, I saved the edits. Now the masked region and the
region outside is shown on the screen. Select the mask region
using the black arrow located on the editor toolbar. Only the
masked region is selected (highlighted in blue). Verify
that the target (on the editor toolbar) is the layer you just created
(i.e., not the masked layer, but the layer which contains both the
masked layer and the area outside). Now copy and paste the masked
region into the outline region. You can verify this by turning
off the masked region layer. Now open the attribute table, there
should
be two rows, one for the masked region, and one for the region outside.
Click on the rows to determine the masked area (i.e.
volcanic area) and change this value to 0. Also choose the row
representing the DEM
boundary and change this value to 1. Stop editing and save the
edits. Unfortunately there is no tool which allows the extraction
outside of a masked area, so this process has to be done instead. Be
sure to search on the current version of ArcGIS to see if this tool is
available for raster datasets, if so this process will not have to be
completed.
Now I have one shapefile with the masked region and the region
outside. Next I wanted to convert this region to a raster. Add the
spatial analyst toolbar and choose convert feature to raster (shown
below).
Here the input feature and the appropriate field is chosen, in this
case the field is Value. Use the cell size of the current DEMs, in this
case that is 10 m and output raster is also specified. Now there is a
raster with 0 and 1 values, representing the masked region and areas
outside respectively. Next the raster needs to be reclassified so the
masked region is represented as NoData and the region outside is 1.
Next, multiply the raster by the difference DEM. To
create the difference DEM, in this case 2001-1954, the
Raster Calculator is used to subtract the datasets. Once this
difference
DEM is created this will be multiplied by the raster we just created (1
and
NoData values). This is done using the Raster Calculator
located on the Spatial Analyst Toolbar. Both tools are shown
below.
Once the above steps are completed we now have a difference DEM with
the volcanic region represented as no data as shown below. The region
shown in step 1 of the overall kriging steps (above) represents the area
outside of the volcanic zone.
Note: When created rasters using reclassify and raster calculator
a layer is normally created called Calculation. Export this file
by right clicking on the layer and choosing Data. Save the file with a
more appropriate file name.
2. Next I
converted each pixel to a point shapefile. But first I changed
the DEM to an integer file instead of a floating point file. Integer
and floating point files are discussed here.
First identity areas on the DEM to determine the number of values to the
right of the decimal point. We need to multiply this value by the DEM
to obtain only integer values. To do this we need the raster calculator
again as shown below.
![Raster Calculator. Type Expression: int ( [test2] * 100 ) Multiply by whatever value is needed to change the values to integers](rastercalc2.JPG)
Here the "int" is the command to change a file to an integer. Here
we took the DEM, in this case test2, and multiplied this by 100 to
change the value to an integer. Next we need to convert each pixel to a
point file. This is done using the Spatial Analyst
Toolbar>Convert>Raster to Feature. Specify the field value
which has elevation values, change the output geometry to point, and
save the output.
Now there should be a point at every pixel location with the elevation
value specified as an integer. Next I added a field in the
attribute table (see step 1 of kriging steps, above). When adding
a field the type has to be specified; in this case the type is float,
which requires
values for precisions and scale. Precision represents the number of
digits stored in a field. The scale is the number of digits to the
right of the decimal place. For this example I chose 6 for
precision, and 2 for scale, but choose accordingly for other
datasets. Once the field is created I chose the field
calculator to calculate values for that field as shown below.
![Right click on Attribute name. Choose Field Calculator. Type the expression: [GRID_CODE]/100](calc_float.jpg)
Once I pressed calculate values, I pressed yes to edit outside of an
edit
session. Using the field calculator I created the Z column to be
Grid_Code/100. The Grid_Code field is the integer value that was
previously created. Now I have a point shapefile which has the correct
elevation values.
3. Next I chose a random sample of the points outside
of the volcanic zone. I needed to interpolate the region within
the volcanic zone and did not want to use all of the points due to an
over representation and processing time. For this study I used a random
sample of 1000. To obtain a random sample the Hawth's tools
extension was downloaded, which is available here.
This toolset currently works with ArcGIS 9.0. If the toolset is
not available for the version you are using, then you will have to find
another method to obtain a random sample, such as Matlab. Using
Hawth's tools a user can choose a random sample and specify the
number of features needed. For this study I created a random sample of 1000 points.
Step 1 in the overall kriging steps (above) shows the 1000 randomly
sampled points overlayed on the difference DEM.
4. 1000 randomly selected points were used in
the kriging analysis to interpolate a surface. Again I highly
suggest researching kriging and geostatistics prior to performing this
analysis, some references can by found here.
I will provide the basic steps, but the values used for the
semivariogram need to be obtained for each data set and an expert
should be consulted to ascertain these. First I will provide a
background for semivariogram and what the definitions mean.
The semivariogram measures the spatial
dependence of points; points which are close together are more alike than
points which are farther apart (Davis, 1986; Maune, 2007). The semivariogram is known only at discrete
points representing distances. To apply the semivariogram for any distance, a
theoretical model needs to be defined which best fits the semivariogram,
therefore creating a continuous function instead of a discrete one (Davies,
1986). A model is chosen to best fit the
semivariogram, which includes measurements for the nugget, sill, and range. The
nugget represents
measurement/independent error and is the deviation from 0 on the y-axis. The sill is the height the semivariogram
reaches when it levels off, located on the y-axis, and the range is the distance where the model first flattens out on the
x-axis. The above parameters are obtained from the semivariogram plot and are
input into the interpolation program using ordinary kriging (OK) to create an interpolated prediction surface. In addition, an error surface was created which analyzes
how accurate the prediction surface is.
Now that we have some basic definitions the fitting of the semivariogram and kriging techniques are shown below:
After choosing the geostatistical wizard the point file needs to be
used for the input data with the appropriate attribute field. There are
various statistical methods which can be used, but kriging is chosen
for this study. In the next menu choose Ordinary kriging, first
the user would create the prediction map and then later would create
the prediction standard error map. The next menu is shown below:
The above menu shows the semivariogram on the left. The semivariogram
is analyzed by a geostatistics expert and the following parameters:
model, range, sill, and nugget are determined and input into the above
menu. When finished with this menu continue through the rest,
unless additional parameters are needed. When finished with the
geostatistical wizard a surface will be created from the interpolation of
the1000 points. The kriging prediction surface and the 1000 randomly
sampled points are shown below:
.
We have now completed step 3 in the overall kriging steps (shown above)
5. Next I clipped the kriging prediction surface to the
masked region, in this case the dome. To recall how to extract an area
by a mask refer here. Once
the kriging prediction surface is clipped to the area of interest the
raster calculator is then used to subtract the kriging surface from the
original DEM difference as shown in the overall kriging steps part 4-6. This
new surface represents the DEM difference with the error removed. For
this study, the kriging prediction surface is the DEM error. More
detailed information is located within the thesis.
6. Block kriging was
used as the final step for the kriging/error analysis.
Prior to block kriging the 1000 points need to be converted to a text
file. The text file should include X, Y, and Z values. To add the X and
Y values I used the XToolsPro toolbar, which is an extension in ArcGIS
9.0. If this extension is not available another method will have
to be used to obtain these values. Exporting a textfile is done by
opening the attribute table, choosing options,
and export. Here you are able to save the file as a text file.
The text file will then be input into the kriging menus.
Block kriging
uses an average number of points within a specified area, which
smoothes the
data. For block kriging the semivariogram calculates the covariance
between
blocks instead of points. To determine the volume change error estimate
we need
to use a comparable density of points. For this study the dome region
has
elevation change outside of the volcanic zone, with a large gap in the
volcanic
zone. In order to interpolate in the volcanic zone, blocks are
needed, which
represent the spacing (area) of the volcanic zone. Otherwise
if 1000 points were used for example, an over representation of
points (or
blocks) would occur outside of the volcanic zone, which would
incorrectly interpolate the
region within the volcanic zone. Therefore blocks were chosen to
represent the
volcanic area. For this study the dome had a block spacing of
2500
m, resulting in an area of 6.25 km2.
Block kriging was performed in Surfer 7.0® by inputting the model
parameters determined from the semivariogram as shown below:
Using Surfer 7.0® the data is gridded using Grid Data available
under the Grid tool. In this menu the x,y, and z
columns are assigned, the kriging method and spacing are also chosen. Using the advanced
options, the model and other kriging parameters are assigned, as well as changing the kriging type to block kriging. . A standard deviation error surface is then created (overall kriging step 7, shown above). These
values were then averaged for the volcanic region, giving a standard error of the estimate. We assume the errors of estimation are
normally distributed; therefore multiplying by 2, we obtain an error estimate
with 95% probability. This value is then
multiplied by the volcanic area to obtain a standard
error for the volume change (overall kriging
step 8, shown above). Again I would suggest consulting a
geostatistics expert before performing the block kriging.
Conclusions
For this study three assumptions were made:
1. That
areas outside of the analysis masks have 0 m elevation change. The region outside of the masked volcanic
region should only be affected by regional uplift, which we assume is
negligible compared to the other DEM errors.
2. DEM error is represented by the deviation
from 0 m elevation change outside of the masked volcanic zone. For example the
average elevation change of 8 m outside of the dome volcanic zone (see thesis) is
assumed to be solely related to DEM error and not to actual elevation change.
3. The errors for the DEM are the same, both inside and
outside, of the volcanic region; thus we can use the elevation changes outside
of the volcanic region to determine an error surface inside the volcanic
region. Since the underlying DEM data was obtained in the same way (aerial
photographs and maps converted to digital contour lines), we can make the
assumption that the same errors should apply to the entire DEM including
volcanic and non-volcanic areas.
Please be aware of these assumptions if you plan on using the same analysis.
Many resources are available on the web related to geostatistics,
kriging, geostatistical analyst, and surfer with some available here.
Some final conclusions and future research ideas are located here.