reinstate.sml

# Reinstate.sml
#
# Script designed to replace interpolated values in a filtered DEM
# Raster with original values if the interpolated value is within a
# tolerance of the original.
#
# Method:
# -------
#
# Requires an unfiltered DEM and a filtered, interpolated DEM as input.
# 1. Calculate difference between unfiltered DEM and filtered DEM
# 2. If difference is less than delta, copy value from unfiltered DEM
#    to output Raster. If its bigger than delta use the value from the
#    filtered DEM.
#
# Terranean Mapping Technologies, Brisbane, Australia: 
# 12/5/2009
#####################################################################
clear();
class RVC_RASTER DEM, FilteredDEM, DEM_reinst;
#Get delta value from user
numeric delta = PopupNum("Threshold vertical distance", 0.4, 0.01, 10, 2);
#Get input rasters
GetInputRasters(DEM,FilteredDEM);
string rvc$ = GetObjectFileName(FilteredDEM);
string name$ = GetObjectName(rvc$,GetObjectNumber(FilteredDEM));
numeric nullValue = NullValue(FilteredDEM);
#Create output DEM
string desc$ = "Reinstate "+FilteredDEM.$Info.Desc;
CreateRaster(DEM_reinst,rvc$,name$+"_edgefill",desc$,NumLins(FilteredDEM),NumCols(FilteredDEM),"32-bit float");
CopySubobjects(FilteredDEM,DEM_reinst, "georef");
#Get output Raster and copy georeference.
GetOutputRaster(DEM_reinst,NumLins(DEM),NumCols(DEM),"32-bit float");
CopySubobjects(DEM,DEM_reinst,"georef");
#Loop through raster and calculate difference between DEM and Filtered
#then set output DEM value accordingly.
foreach FilteredDEM{
if (DEM - FilteredDEM < delta) DEM_reinst = DEM;
else DEM_reinst = FilteredDEM;
}
#Clean up
CreateHistogram(DEM_reinst);
CreatePyramid(DEM_reinst);
#EOF