######################################################################################################
#
# SRTMfill.sml ("Stand Alone" Script)
#
######################################################################################################
#
# Authors:
# Dave Breitwisch, MicroImages, Inc
#
# Started: 07 April 2005
# Updated: 28 March 2006
#
# Requires: TNTmips2005:71 or later
#
numeric version = round(_context.Version * 10);
if (version < 71) {
PopupMessage("This version requires that you are using TNTmips2005:71 or later");
Exit();
}
#
######################################################################################################
#
# Description of Script:
#
# This is a "Stand Alone" script whose purpose is to fill in null cell values in the released data
# from the Shuttle Radar Topography Mission flown in February of 2000.
# This current version of the script:
# fills single cells by interpolation,
# fills null voids that are bounded by a constant value, and
# fills voids using elevation data from reference DEM rasters provided by the user.
#
#
# Void Filling Procedures:
#
# Single Cell Interpolation - The user can specify the number of surrounding cells that are
# required (minimum = 6) to be present for the single cell to be filled. If the single null
# cell passes this test, the average of all the surrounding non-null cells is used. This procedure
# is called multiple times until there are no remaining single null cells that pass the test.
#
# Constant Boundary Filling - The algorithm used for this procedure searches for null voids that
# have a constant value around their edges. When a null void is found, the algorithm
# traverses around the boundary of the void. If all the surrounding values are the same,
# the algorithm fills in it's traversed path with that constant value. The remaining cells
# inside of this path are filled in when the algorithm detects this new null void. This
# is especially useful for filling in voids within a water surface.
# NOTE: If a null void covers a real area of different elevation, but the surrouding cells
# are the same value, the void is filled in with the constant value. This is the
# case when an SRTM void includes an island within a horizontal water body.
#
# Reference DEM Filling - This allows the user to use multiple reference DEMs, assumed to be
# accurate and properly georeferenced, to fill the null voids. The DEM data is assumed to be
# in meters unless the proper cell value scale is set for the raster. Additional options such
# as a maximum elevation difference and blending are available.
#
# Maximum Elevation Difference - For each reference DEM, the values of the SRTM cells are
# compared to the matching georeferenced values of the DEM. If the differnece between the
# values is greater then specified by the user, the SRTM cell is replaced by null and will
# eventually be filled in with the reference DEM value. This is done under the assumption
# that the reference DEMs are more accurate.
#
# Blending - This procedure uses the value specified by the user to blend the areas around
# the null voids filled in by the reference DEMs. The value represents the distance away
# from the boundary of the null void (in number of cells) over which to blend. Blending
# is performed using nested 1-cell-wide buffer zones. The final SRTM cell values in each
# buffer zone are linear weighted averages of the original SRTM cell value and the reference
# DEM cell value with the weighting factors determined by the number of buffer zones and
# the position of the current buffer zone in the nesting sequence.
#
# Example: blending over 4 nested buffer zones around each void
# Zone 1 : cell = SRTM value + 4/5 * diffence between the values (Boundary of void)
# Zone 2 : cell = SRTM value + 3/5 * diffence between the values
# Zone 3 : cell = SRTM value + 2/5 * diffence between the values
# Zone 4 : cell = SRTM value + 1/5 * diffence between the values (Furthest from void)
#
# An additional option is given to the user to save the raster and vector objects created
# by this filling process. These can be used for data verification. These objects are made
# after the Single Cell Interpolation and Constant Boundry Filling operations are done, if
# chosen by the user. The following objects are created:
# BinNullVoids (Binary Raster): 1 = Matching cell to be filled in by reference DEM value
# BufferR# (Binary Raster): 1 = Matching cell is in Zone number #
# NullVoids (Vector): Shows boundaries of BinNullVoids Raster
# BufferV# (Vector): Shows outer boundaries of the matching BufferR# Raster
#
######################################################################################################
######################################################################################################
# Global Variable Declarations
######################################################################################################
class XMLDOC mainDOC, progressDOC;
class XMLNODE mainNODE, progressNODE;
class GUI_DLG mainDLG, progressDLG;
numeric errXML, errDLG;
string mainxmlCODE$, progressxmlCODE$;
string InputRasterFilename$, InputRasterObjectName$;
string OutputRasterFilename$, OutputRasterObjectName$;
string LogFilename$;
class STRINGLIST RefDEMFilenames$, RefDEMObjectNames$;
class RASTER inR;
class RASTER outR;
numeric rows, cols;
string initType;
numeric numValidNulls;
class RASTER RefR;
numeric RefDataScale;
class RASTER binVoidR;
class VECTOR nullVoidsV;
array numeric element[1,1];
array numeric buffinodes[1];
###########################
# User Selected Values
###########################
numeric cellBufferDist;
numeric threshold;
numeric numRequiredAround;
numeric numReferenceDEMs;
string TempSRTMfile;
######################################################################################################
######################################################################################################
# Dialog specification in XML
######################################################################################################
mainxmlCODE$ = '
';
######################################################################################################
######################################################################################################
# Dialog specification in XML
######################################################################################################
progressxmlCODE$ = '
';
######################################################################################################
######################################################################################################
# Start of the Script
######################################################################################################
clear();
proc OnCreateDialog(); #Function Prototype
OnCreateDialog();
######################################################################################################
######################################################################################################
# Script Functions and Procedures
######################################################################################################
######################################################################################################
# Returns the number of NULL cells that are in the raster
func NumNullCellsInRaster(class RASTER rast) {
local numeric i, j, nullCount = 0;
for i = 1 to NumLins(rast) {
for j = 1 to NumCols(rast) {
if (IsNull(rast[i,j])) {
nullCount++;
}
}
}
return nullCount;
}
######################################################################################################
# Return the number of cells surrounding the location of the cell passed that are null.
func NumNullsAround(numeric i, numeric j) {
local numeric count = 0;
# Add 1 to count if specified neighboring cell is a valid Raster cell and is Null
if ((i != 1) && (j != cols) && IsNull(outR[i-1,j+1])) count++;
if ((i != 1) && IsNull(outR[i-1,j ])) count++;
if ((i != 1) && (j != 1) && IsNull(outR[i-1,j-1])) count++;
if ( (j != cols) && IsNull(outR[i ,j+1])) count++;
if ( (j != 1) && IsNull(outR[i ,j-1])) count++;
if ((i != rows) && (j != cols) && IsNull(outR[i+1,j+1])) count++;
if ((i != rows) && IsNull(outR[i+1,j ])) count++;
if ((i != rows) && (j != 1) && IsNull(outR[i+1,j-1])) count++;
return count;
}
######################################################################################################
# Returns the average value of all surrounding cells that are not null
func InterpCellValue(numeric i, numeric j) {
local numeric sumvalue = 0;
local numeric numCells = 0;
# Check if the neighboring cell is valid and not null. If ok, add to sum and increment the count
if ((i != 1) && (j != cols) && (!IsNull(outR[i-1,j+1]))) {
sumvalue += outR[i-1,j+1];
numCells++;
}
if ((i != 1) && (!IsNull(outR[i-1,j]))) {
sumvalue += outR[i-1,j];
numCells++;
}
if ((i != 1) && (j != 1) && (!IsNull(outR[i-1,j-1]))) {
sumvalue += outR[i-1,j-1];
numCells++;
}
if ((j != cols) && (!IsNull(outR[i,j+1]))) {
sumvalue += outR[i,j+1];
numCells++;
}
if ((j != 1) && (!IsNull(outR[i,j-1]))) {
sumvalue += outR[i,j-1];
numCells++;
}
if ((i != rows) && (j != cols) && (!IsNull(outR[i+1,j+1]))) {
sumvalue += outR[i+1,j+1];
numCells++;
}
if ((i != rows) && (!IsNull(outR[i+1,j]))) {
sumvalue += outR[i+1,j];
numCells++;
}
if ((i != rows) && (j != 1) && (!IsNull(outR[i+1,j-1]))) {
sumvalue += outR[i+1,j-1];
numCells++;
}
# Return the average of the sum
return sumvalue/numCells;
}
######################################################################################################
# This process traverses through the SRTM raster and looks for single null cells that have less then
# the allowed number of nulls surrounding them. For these cells the function interpolates a value
# based on the average of the surrounding real cell values.
func CheckNegligibleNulls( numeric allowedNulls ) {
local numeric cellsChanged = 0;
SetStatusMessage("Interpolating Negligible Null Cells");
local numeric i, j;
# Traverse through the raster
for i = 1 to rows {
SetStatusBar(i, rows);
for j = 1 to cols {
# Check if the cell is null
if (IsNull(outR[i, j])) {
# Check if it has less then the number of surrounding null cells allowed
if ( NumNullsAround(i, j) <= allowedNulls ) {
# Interpolate the cell value and update the report counter
outR[i, j] = InterpCellValue(i, j);
cellsChanged++;
}
}
}
}
print(" Number of single null cells changed in this pass: " + NumToStr(cellsChanged) + "\n");
return cellsChanged;
}
######################################################################################################
# This procedure creates a binary raster that symbolizes the overlap of the input SRTM raster and the
# reference DEM. The value of 1 indicates that the SRTM has a null cell. A vector is also created
# that contains the boundaries of the null cells (value of 1) in the binary raster.
proc CreateNullRasterVector( ) {
print("Creating Void Mask and Void Boundries\n");
CreateRaster(binVoidR, TempSRTMfile, "BinNullVoids", "Binary Raster showing Null Voids", rows, cols, "binary");
CopySubobjects(outR, binVoidR, "GEOREF");
SetNull(binVoidR, 0);
# Set the cell to 1 if the matching SRTM raster cell is null, else set it to 0.
local numeric i, j, I = 1, J = 1;
for i = 1 to rows {
for j = 1 to cols {
# Offset the line and column counters to match the two rasters
if ( IsNull(outR[i, j]) ) {
binVoidR[I,J] = 1;
} else {
binVoidR[I,J] = 0;
}
J++;
}
I++;
J = 1;
}
# Create the vector with boundaries around the 1's in the binary raster
CreateVector(nullVoidsV, TempSRTMfile, "NullVoids", "Boundary of Null Voids (after single cell pass)");
nullVoidsV = RasterToVectorBound(binVoidR, "DoImpliedGeoref");
# Save the object numbers so they can be reopened later
ResizeArrayClear(buffinodes, (cellBufferDist+1)*2);
buffinodes[1] = GetObjectNumber(binVoidR);
buffinodes[2] = GetObjectNumber(nullVoidsV);
CloseRaster(binVoidR);
CloseVector(nullVoidsV);
}
######################################################################################################
# This procedure traverses through the SRTM raster and replaces any null values with the corresponding
# cell value from the reference DEM. The SRTM cell will remain null if the current reference DEM does
# not overlap or if the corresponding reference DEM cell is null.
proc FillWithRefRaster() {
print("Filling Voids with Reference Data\n");
local class GEOREF outG, refG;
local class POINT2D point;
outG = GetLastUsedGeorefObject(outR);
refG = GetLastUsedGeorefObject(RefR);
# create a transparm if the CRS's are different
local numeric same = 1;
if (outG.CoordRefSys.Name != refG.CoordRefSys.Name) {
same = 0;
local class TRANSPARM trans;
trans.InputCoordRefSys = outG.CoordRefSys;
trans.OutputCoordRefSys = refG.CoordRefSys;
}
SetStatusMessage("Filling holes with Reference Data");
# traverse through the overlapping area
local numeric i, j;
for i = 1 to rows {
SetStatusBar(i, rows);
for j = 1 to cols {
# If the cell is null, find the matching cell in the reference DEM
if ( IsNull(outR[i,j]) ) {
point = ObjectToMap(outR, j, i, outG);
# Convert to the correct CRS if they are not the same
if (!same) {
point = trans.ConvertPoint2DFwd(point);
}
# convert map coordinates to SRTM raster line/column
point = MapToObject(refG, point.x, point.y, RefR);
point.x = round(point.x);
point.y = round(point.y);
# Check to make sure the point is valid, then replace it
if ( point.x >= 1 and point.x <= RefR.$Info.NumCols) {
if ( (point.y >= 1) && (point.y <= RefR.$Info.NumLins) ) {
outR[i,j] = RefR[point.y, point.x] * RefDataScale;
}
}
}
}
}
}
######################################################################################################
# This procedure blends the SRTM data with the Reference data that was used to fill the null cells
# based on the distance from the null voids that the user specifies.
proc BlendBufferZone( class RASTER zoneBufferR, numeric zoneNum ) {
print(" Buffering in Zone " + NumToStr(zoneNum) + " of " + NumToStr(cellBufferDist) + "\n");
local class POINT2D point;
local class GEOREF outG, refG;
local numeric I = 1, J = 1;
outG = GetLastUsedGeorefObject(outR);
refG = GetLastUsedGeorefObject(RefR);
# create a transparm if the CRS's are different
local numeric same = 1;
if (outG.CoordRefSys.Name != refG.CoordRefSys.Name) {
same = 0;
local class TRANSPARM trans;
trans.InputCoordRefSys = outG.CoordRefSys;
trans.OutputCoordRefSys = refG.CoordRefSys;
}
SetStatusMessage("Computing new values in Buffer Zone");
# traverse through the overlapping area
local numeric m, n, buffFactor;
for m = 1 to rows {
SetStatusBar(m, rows);
for n = 1 to cols {
# Start the blending process if the cell is in the buffer zone
if ( zoneBufferR[I, J] == 1 ) {
point = ObjectToMap(outR, n, m, outG);
# Convert point to the correct CRS if not the same
if (!same) {
point = trans.ConvertPoint2DFwd(point);
}
point = MapToObject(refG, point.x, point.y, RefR);
point.x = round(point.x);
point.y = round(point.y);
# Make sure the matching cell in the DEM is valid
if ( point.x >= 1 and point.x <= RefR.$Info.NumCols) {
if ( (point.y >= 1) && (point.y <= RefR.$Info.NumLins) ) {
# Assign the output cell the value of the SRTM plus the linear differnce based
# on what zone the cells are in.
buffFactor = (cellBufferDist - zoneNum + 1) * 1.0 / (cellBufferDist * 1.0 + 1);
outR[m,n] = outR[m,n] + (buffFactor) * ((RefR[point.y, point.x] * RefDataScale) - outR[m,n]);
}
}
}
J++;
}
I++;
J = 1;
}
}
######################################################################################################
# This procedure uses the binVoidR raster created early along with the morphological functions to grow
# the voids larger based on the zone number. A boundary vector is created with each new raster.
proc CreateVectRastBufferZoneInit(numeric zone) {
SetStatusMessage("Computing initial void Vector and Raster");
local string zone$ = NumToStr(zone);
print(" Growing Buffer Raster " + zone$ + " of " + NumToStr(cellBufferDist) + "\n");
# Create a raster the size of binVoidR and dilate the null void based on the zone value
local class RASTER buffR;
local numeric tempNumLins = NumLins(binVoidR);
local numeric tempNumCols = NumCols(binVoidR);
CreateRaster(buffR, TempSRTMfile, "BufferR"+zone$, "Buffer Zone "+zone$, tempNumLins, tempNumCols, "binary");
CopySubobjects(binVoidR, buffR, "GEOREF");
SetNull(buffR, 0);
MorphDilation(binVoidR, buffR, element, (zone*2)+1, (zone*2)+1);
# Create a vector that bounds each of the null voids
local class VECTOR buffV;
CreateVector(buffV, TempSRTMfile, "BufferV"+zone$, "Boundary of Buffer Zone "+zone$);
buffV = RasterToVectorBound(buffR, "DoImpliedGeoref");
# Save the inodes for future use
buffinodes[(zone*2)+1] = GetObjectNumber(buffR);
buffinodes[(zone*2)+2] = GetObjectNumber(buffV);
CloseVector(buffV);
CloseRaster(buffR);
}
######################################################################################################
# This procedure compares each of the binary rasters with the raster that is "grown 1 size smaller" so
# the resulting raster will only contain the zone boundaries. All else will be 0's.
proc ComputeBufferZoneRasters( numeric zone ) {
print(" Creating Buffer Zone Raster " + NumToStr(zone) + " of " + NumToStr(cellBufferDist) + "\n");
# Open the working raster and the "1 less grown" raster
local string name;
local class RASTER grownR;
local class RASTER prevR;
name = GetObjectName(TempSRTMfile, buffinodes[(zone*2)+1]);
OpenRaster(grownR, TempSRTMfile, name);
name = GetObjectName(TempSRTMfile, buffinodes[((zone-1)*2)+1]);
OpenRaster(prevR, TempSRTMfile, name);
local numeric tempNumLins = NumLins(grownR);
local numeric tempNumCols = NumCols(grownR);
SetStatusMessage("Computing Buffer Zone Raster");
# Traverse through the working raster
local numeric r, c;
for r = 1 to tempNumLins {#rows {
SetStatusBar(r, tempNumLins);
for c = 1 to tempNumCols {#cols {
# Do and XOR with the two rasters
if ( (IsNull(prevR[r,c])) && (grownR[r,c] == 1) ) {
grownR[r,c] = 1;
} else {
grownR[r,c] = 0;
}
}
}
CloseRaster(grownR);
if ( zone != 1 ) { # Do not close binVoidR
CloseRaster(prevR);
}
}
######################################################################################################
# This procedure controls the process of blending the reference raster data that was used to fill the
# null voids in the SRTM data.
proc BlendWithRefRaster() {
# Create an element large enough to satisfy the users number of blending zones. This is used in
# the morphological dilation procedures
ResizeArrayClear(element, (cellBufferDist*2)+1, (cellBufferDist*2)+1);
numeric k, i;
for k = 1 to (cellBufferDist*2)+1 {
for i = 1 to (cellBufferDist*2)+1 {
element[k,i] = 1;
}
}
# Create the initial binary rastes and grow them according the the zone value
print("Using Raster Dilation to produce buffer Rasters \n");
for k = 1 to cellBufferDist {
CreateVectRastBufferZoneInit(k);
}
# Do an XOR with each of the raster with the "1 lesser grown" raster. Produces binary zones.
print("Using buffer Rasters to create Buffer Zones \n");
for k = cellBufferDist to 1 step -1 {
ComputeBufferZoneRasters(k);
}
print("Buffering around Voids to match data\n");
# Call the procedure that will blend the reference DEM and SRTM data using the zone rasters
local class RASTER blendR;
local string name;
for k = 1 to cellBufferDist {
name = GetObjectName(TempSRTMfile, buffinodes[(k*2)+1]);
OpenRaster(blendR, TempSRTMfile, name);
BlendBufferZone(blendR, k);
CloseRaster(blendR);
}
}
######################################################################################################
# This procedure scans through the overlapping SRTM and reference DEM rasters and removes the SRTM
# value if the elevation difference is greater then the threshold specified by the user. It is then
# replaced with the reference DEM data in a later procedure.
proc RemoveThresholdViolations( ) {
local class POINT2D point;
local class GEOREF outG, refG;
print("Removing Threshold Violations to Referance Raster \n");
outG = GetLastUsedGeorefObject(outR);
refG = GetLastUsedGeorefObject(RefR);
# Create a transparm if the CRSs are not the same
local numeric same = 1;
if (outG.CoordRefSys.Name != refG.CoordRefSys.Name) {
same = 0;
local class TRANSPARM trans;
trans.InputCoordRefSys = outG.CoordRefSys;
trans.OutputCoordRefSys = refG.CoordRefSys;
}
SetStatusMessage("Removing values that violate the specified Threshold to the reference data\n");
# Traverese the SRTM raster
local numeric numViolations = 0;
local numeric m, n;
for m = 1 to rows {
SetStatusBar(m, rows);
for n = 1 to cols {
point = ObjectToMap(outR, n, m, outG);
# Convert the point to the SRTM CRS if not the same
if (!same) {
point = trans.ConvertPoint2DFwd(point);
}
point = MapToObject(refG, point.x, point.y, RefR);
point.x = round(point.x);
point.y = round(point.y);
# Check if the point is valid
if ( point.x >= 1 and point.x <= RefR.$Info.NumCols) {
if ( (point.y >= 1) && (point.y <= RefR.$Info.NumLins) ) {
# Check if the difference in elevation values is greater then the threshold
if ( abs(outR[m,n] - (RefR[point.y, point.x] * RefDataScale)) > threshold) {
# Set the cell value to null and update the count
if (HasNullValue(outR)) {
outR[m,n] = NullValue(outR);
}
if (HasNullMask(outR)) {
outR[m,n] = null;
}
numViolations++;
}
}
}
}
}
print(" Total violations replaced with null: " + NumToStr(numViolations) + "\n");
}
######################################################################################################
# Is true if all the surrouding cells from the specified cell are null or equal to the passed value.
func AreAllSurroundingSameAsValueOrNull(numeric i, numeric j, numeric value) {
local numeric valid = true;
# Checks if cell exists. Returns false if the existing cell is Null
if ((outR[i-1,j+1] != value) & !IsNull(inR[i-1,j+1])) valid = false;
if ((outR[i-1,j ] != value) & !IsNull(inR[i-1,j ])) valid = false;
if ((outR[i-1,j-1] != value) & !IsNull(inR[i-1,j-1])) valid = false;
if ((outR[i ,j+1] != value) & !IsNull(inR[i ,j+1])) valid = false;
if ((outR[i ,j-1] != value) & !IsNull(inR[i ,j-1])) valid = false;
if ((outR[i+1,j+1] != value) & !IsNull(inR[i+1,j+1])) valid = false;
if ((outR[i+1,j ] != value) & !IsNull(inR[i+1,j ])) valid = false;
if ((outR[i+1,j-1] != value) & !IsNull(inR[i+1,j-1])) valid = false;
return valid;
}
######################################################################################################
# Is true if all the surrounding cells from the specified cell are equal to the passed value.
func AreAllSurroundingSameAsValue(numeric i, numeric j, numeric value) {
# Check each surrouding cell
if ( IsNull(outR[i-1,j-1]) ) {
return false;
} else if ( outR[i-1,j-1] != value ) {
return false;
}
if ( IsNull(outR[i-1,j ]) ) {
return false;
} else if ( outR[i-1,j ] != value ) {
return false;
}
if ( IsNull(outR[i-1,j+1]) ) {
return false;
} else if ( outR[i-1,j+1] != value ) {
return false;
}
if ( IsNull(outR[i ,j-1]) ) {
return false;
} else if ( outR[i ,j-1] != value ) {
return false;
}
if ( IsNull(outR[i ,j+1]) ) {
return false;
} else if ( outR[i ,j+1] != value ) {
return false;
}
if ( IsNull(outR[i+1,j-1]) ) {
return false;
} else if ( outR[i+1,j-1] != value ) {
return false;
}
if ( IsNull(outR[i+1,j ]) ) {
return false;
} else if ( outR[i+1,j ] != value ) {
return false;
}
if ( IsNull(outR[i+1,j+1]) ) {
return false;
} else if ( outR[i+1,j+1] != value ) {
return false;
}
return true;
}
######################################################################################################
# This is a recursive function. It is called when a null cell is found. It uses the sequence of
# travel specified below when traversing around the raster. As it travels around the boundry of the
# null void, it checks to see if the surrounding cells are the same constant value as the cell value
# to the left of the initial null cell. If the function ends up at the same point it started and did
# not fail any of it's surrounding cell checks, it goes back through the same path that it traveled
# and replaces each cell with the constant value. If a void still exists within this traversed path,
# it will be filled in as the original function that calls this recursive function traverses the SRTM
# raster.
func RecursiveBorderCheck(numeric i, numeric j, numeric start_i, start_j, string traveled, numeric value, numeric firstCall) {
# If the current position one cell to the right of starting cell and it is not the first call of
# the function, it means that it has gone all the way around the border and each border cell has
# the same constant value. Returns true so it recusively fills in each of the cells that the
# function was called on with the constant value.
if (( i == start_i) & (j-1 == start_j) & (!firstCall)) {
return true;
}
# Check if the surrounding cells pass the constant value or null test
if (AreAllSurroundingSameAsValueOrNull(i, j, value)) {
# Determine which direction to try and travel based on the sequence below:
#
# Sequence of travel:
# If the function was called by traveling to the:
# RIGHT - try in order of - UP, RIGHT, DOWN, LEFT
# DOWN - try in order of - RIGHT, DOWN, LEFT, UP
# LEFT - try in order of - DOWN, LEFT, UP, RIGHT
# UP - try in order of - LEFT, UP, RIGHT, DOWN
# Check which direction the calling funtion traveled to the current cell
if (traveled == "RIGHT") {
# If the cell UP is null, go that way
if (IsNull(outR[i-1, j])) {
# Call this function and go that way
if (RecursiveBorderCheck(i-1, j, start_i, start_j, "UP", value, 0)) {
# If it returned true, it is traversing back after passing the constant border test
outR[i,j] = value;
return true;
} else {
# Failed the constant border test, just exit out of the function
return false;
}
# If the cell RIGHT is null, go that way
} else if (IsNull(outR[i,j+1])) {
if (RecursiveBorderCheck(i, j+1, start_i, start_j, "RIGHT", value, 0)) {
outR[i,j] = value;
return true;
} else {
return false;
}
# If the cell DOWN is null, go that way
} else if (IsNull(outR[i+1,j])) {
if (RecursiveBorderCheck(i+1, j, start_i, start_j, "DOWN", value, 0)) {
outR[i,j] = value;
return true;
} else {
return false;
}
# If the cell LEFT is null, go that way
} else if (IsNull(outR[i,j-1])) {
if (RecursiveBorderCheck(i, j-1, start_i, start_j, "LEFT", value, 0)) {
outR[i,j] = value;
return true;
} else {
return false;
}
}
}
# Comments similar to the above
if (traveled == "DOWN") {
if (IsNull(outR[i, j+1])) {
if (RecursiveBorderCheck(i, j+1, start_i, start_j, "RIGHT", value, 0)) {
outR[i,j] = value;
return true;
} else {
return false;
}
} else if (IsNull(outR[i+1,j])) {
if (RecursiveBorderCheck(i+1, j, start_i, start_j, "DOWN", value, 0)) {
outR[i,j] = value;
return true;
} else {
return false;
}
} else if (IsNull(outR[i,j-1])) {
if (RecursiveBorderCheck(i, j-1, start_i, start_j, "LEFT", value, 0)) {
outR[i,j] = value;
return true;
} else {
return false;
}
} else if (IsNull(outR[i-1,j])) {
if (RecursiveBorderCheck(i-1, j, start_i, start_j, "UP", value, 0)) {
outR[i,j] = value;
return true;
} else {
return false;
}
}
}
# Comments similar to the above
if (traveled == "LEFT") {
if (IsNull(outR[i+1, j])) {
if (RecursiveBorderCheck(i+1, j, start_i, start_j, "DOWN", value, 0)) {
outR[i,j] = value;
return true;
} else {
return false;
}
} else if (IsNull(outR[i,j-1])) {
if (RecursiveBorderCheck(i, j-1, start_i, start_j, "LEFT", value, 0)) {
outR[i,j] = value;
return true;
} else {
return false;
}
} else if (IsNull(outR[i-1,j])) {
if (RecursiveBorderCheck(i-1, j, start_i, start_j, "UP", value, 0)) {
outR[i,j] = value;
return true;
} else {
return false;
}
} else if (IsNull(outR[i,j+1])) {
if (RecursiveBorderCheck(i, j+1, start_i, start_j, "RIGHT", value, 0)) {
outR[i,j] = value;
return true;
} else {
return false;
}
}
}
# Comments similar to the above
if (traveled == "UP") {
if (IsNull(outR[i, j-1])) {
if (RecursiveBorderCheck(i, j-1, start_i, start_j, "LEFT", value, 0)) {
outR[i,j] = value;
return true;
} else {
return false;
}
} else if (IsNull(outR[i-1,j])) {
if (RecursiveBorderCheck(i-1, j, start_i, start_j, "UP", value, 0)) {
outR[i,j] = value;
return true;
} else {
return false;
}
} else if (IsNull(outR[i,j+1])) {
if (RecursiveBorderCheck(i, j+1, start_i, start_j, "RIGHT", value, 0)) {
outR[i,j] = value;
return true;
} else {
return false;
}
} else if (IsNull(outR[i+1,j])) {
if (RecursiveBorderCheck(i+1, j, start_i, start_j, "DOWN", value, 0)) {
outR[i,j] = value;
return true;
} else {
return false;
}
}
}
} else {
# Failed the constant surrounding values check, exit out.
return false;
}
}
######################################################################################################
# This procedure traverses through the SRTM raster, and if it finds a null value, it calls a recursive
# function to travel around the border and fill in the void if all surrouding cells have the same
# cell value.
proc FillVoidsWithConstantBoundaryValues( ) {
local numeric i, j;
local numeric value, total = 0;
print("Filling null voids that have a constant value on all borders\n");
SetStatusMessage("Filling null voids with constant boundaries");
# Fill any voids that have a constant boundary surrouding. Will "errode" in one pixel.
# Continuing passes through the raster will keep filling in one pixel at a time. This
# may leave a single cell remaining null based on original shape of the void.
# WARNING: This may partialy or completely cover an island with water if the bounding
# value around the island (void) is constant.
# Traverse through the SRTM raster, except for the outside edge. The algorighm can't determine if
# the "real world elevation" value beyond the raster has the same value since the data does not
# exist in this raster.
for i = 2 to rows-1 {
SetStatusBar(i, rows-1);
for j = 2 to cols-1 {
if (IsNull(outR[i,j])) {
# If the cell to the left and up are not null, get the value and start the processes
if ((!IsNull(outR[i,j-1])) & (!IsNull(outR[i-1,j]))) {
value = outR[i,j-1];
# Do an initial check of the surrouding cells and call the function if it passes
if (AreAllSurroundingSameAsValueOrNull(i, j, value)) {
# Keep a report for each successful filling. (multiple sucesses may be recorded
# since the void is filled in one boudary ring at a time)
total += RecursiveBorderCheck(i, j, i, j-1, "RIGHT", value, 1);
}
}
}
}
}
print(" Number of null voids filled: " + NumToStr(total) + "\n");
total = 0;
# Fill any single cells remaining that have constant values surrounding it. This is done to fill
# in any single cells that the previous recursive function missed since it cannot return true on
# the first call to itself.
for i = 1 to rows-1 {
SetStatusBar(i, rows-1);
for j = 1 to cols-1 {
# Check if the cell is null and the cell to the left is not null
if ( (IsNull(outR[i,j])) & (!IsNull(outR[i,j-1])) ) {
value = outR[i,j-1];
# Check if all surrounding cells have the same value
if (AreAllSurroundingSameAsValue(i, j, value)) {
# Fill in the cell and update the count
outR[i,j] = value;
total++;
}
}
}
}
print(" Number of single null cells filled: " + NumToStr(total) + "\n");
}
######################################################################################################
# Assign the global variables the values set by the user in the dialog
proc SetUserSpecifiedGlobalValues( ) {
numRequiredAround = mainDLG.GetCtrlByID("id_NUM_REQUIRED_SURROUNDING").GetValueNum();
numReferenceDEMs = mainDLG.GetCtrlByID("id_NUM_REF_DEM").GetValueNum();
cellBufferDist = mainDLG.GetCtrlByID("id_NUM_BUFFER_ZONE").GetValueNum();
threshold = mainDLG.GetCtrlByID("id_THRESHOLD_VALUE").GetValueNum();
}
######################################################################################################
# This is the main procedure that is called when the user presses RUN with the appropriate settings
# set in the main dialog. It calls functions based on the processes that the user has requested be
# run on the data provided.
proc RunVoidFillingProcess( ) {
local numeric nullCellsExist = 1;
local class GUI_CTRL_TOGGLEBUTTON saveToggle = mainDLG.GetCtrlByID("id_SAVE_TEMP_FILE_TOGGLE");
SetUserSpecifiedGlobalValues();
print("Starting...\n");
# Check input SRTM raster to see if it has any null cells
OpenRaster(inR, InputRasterFilename$, InputRasterObjectName$);
if (!HasNull(inR)) {
PopupMessage("This raster does not have any recognized Null values");
print("Done!\n");
return;
}
# If the temporary working project file is not specified by the user, create one
if (!saveToggle.GetValue()) {
TempSRTMfile = _context.ScriptDir + "/SRTMvoidfillingTempFile.rvc";
if (fexists(TempSRTMfile)) {
DeleteFile(TempSRTMfile);
}
CreateProjectFile(TempSRTMfile, "Temporary working directory for SRTM void filling SML script");
}
OpenRaster(outR, OutputRasterFilename$, OutputRasterObjectName$);
# Make a copy of the input raster. All procedures will only modify this copy
DeleteSubobjects(outR);
CopySubobjects(inR, outR, "ALL");
if(HasNullValue(inR)) {
SetNull(outR, NullValue(inR));
}
outR = inR;
# Flush the raster object to recognize the null mask if it has one
CloseRaster(outR);
OpenRaster(outR, OutputRasterFilename$, OutputRasterObjectName$);
# Call the procedure to fill single null cells if selected
if (mainDLG.GetCtrlByID("id_FILL_SINGLE_NULLS_TOGGLE").GetValueNum()) {
local numeric negNullTotal = 0;
local numeric temp = 0;
print("Filling single null cells\n");
# Traverse the SRTM raster until no changes are made
while ( (temp = CheckNegligibleNulls(8 - numRequiredAround )) > 0 ) {
negNullTotal = negNullTotal + temp;
}
print(" Total number changed: " + NumToStr(negNullTotal) + "\n");
# Check to see if there are any remaining null cells in the raster
if (NumNullCellsInRaster(outR) == 0) {
print("There are no remaining null cells in the raster. Remaining filling methods will be skipped.\n");
nullCellsExist = 0;
}
}
# Call the procedure to fill null voids with constant boundaries if selected
if (mainDLG.GetCtrlByID("id_FILL_CONSTANT_NULL_VOIDS_TOGGLE").GetValueNum() and nullCellsExist) {
FillVoidsWithConstantBoundaryValues();
# Check to see if there are any remaining null cells in the raster
if (NumNullCellsInRaster(outR) == 0) {
print("There are no remaining null cells in the raster. Remaining filling methods will be skipped.\n");
nullCellsExist = 0;
}
}
# Call the procedure to fill null voids with additional data if selected
if (mainDLG.GetCtrlByID("id_FILL_VOIDS_WITH_REF_DEM_TOGGLE").GetValueNum() and nullCellsExist) {
local numeric i;
# For each reference raster that was selected
for i = 0 to numReferenceDEMs-1 {
if (nullCellsExist) {
OpenRaster(RefR, RefDEMFilenames$[i], RefDEMObjectNames$[i]);
RefDataScale = RefR.$Info.DataScale;
# Call the procedure to replace cells that violate the maximum elevation difference
if (mainDLG.GetCtrlByID("id_USE_THRESHOLD_TOGGLE").GetValueNum()) {
RemoveThresholdViolations();
}
# Create working raster and vectors and fill holes with reference data
CreateNullRasterVector();
FillWithRefRaster();
# Call the procedure to blend around the null voids if selected
if (mainDLG.GetCtrlByID("id_USE_BUFFER_ZONE_TOGGLE").GetValueNum()) {
BlendWithRefRaster();
}
CloseRaster(RefR);
# Check to see if there are any remaining null cells in the raster
if (NumNullCellsInRaster(outR) == 0) {
print("There are no remaining null cells in the raster. Remaining filling methods will be skipped.\n");
nullCellsExist = 0;
}
}
}
}
# Update subobjects for new SRTM raster
CreateHistogram(outR);
CreatePyramid(outR);
print("Void-filled SRTM Raster saved as: " + outR.$Info.Name);
print("Located in project file: " + outR.$Info.Filename + "\n");
CloseRaster(inR);
CloseRaster(outR);
# Delete the working directory if not selected to save
if (!saveToggle.GetValue()) {
DeleteFile(TempSRTMfile);
}
print("Done!\n");
return;
}
######################################################################################################
######################################################################################################
# Dialog Related Functions
######################################################################################################
######################################################################################################
# Called when status dialog is created to start the void filling processes
proc OnProgressOpen( ) {
RunVoidFillingProcess();
progressDLG.Close(0);
}
######################################################################################################
# Called when status dialog is created to check if the selected input has any valid null cells
proc OnVerifyOpen( ) {
# Check to see if there are any null cells in the raster
local class RASTER checkR;
OpenRaster(checkR, InputRasterFilename$, InputRasterObjectName$);
numValidNulls = NumNullCellsInRaster(checkR);
CloseRaster(checkR);
progressDLG.Close(0);
}
######################################################################################################
# Called when "Run..." is pressed. Checks to make sure settings are valid before starting the filling
# procedures.
proc OnRunPressed( ) {
local class GUI_CTRL_PUSHBUTTON buttonRun = mainDLG.GetCtrlByID("id_RUN_BUTTON");
buttonRun.SetEnabled(0);
local class GUI_CTRL_TOGGLEBUTTON toggleSingle = mainDLG.GetCtrlByID("id_FILL_SINGLE_NULLS_TOGGLE");
local class GUI_CTRL_TOGGLEBUTTON toggleConst = mainDLG.GetCtrlByID("id_FILL_CONSTANT_NULL_VOIDS_TOGGLE");
local class GUI_CTRL_TOGGLEBUTTON toggleDEM = mainDLG.GetCtrlByID("id_FILL_VOIDS_WITH_REF_DEM_TOGGLE");
if (!toggleSingle.GetValue() and !toggleConst.GetValue() and !toggleDEM.GetValue()) {
string msg = "Error: Please select at least one void filling procedure: \n\n";
msg = msg + " 1) Interpolate single null cells\n";
msg = msg + " 2) Fill null voids surrounded by constant value\n";
msg = msg + " 3) Fill null voids using reference DEM rasters";
PopupMessage(msg);
buttonRun.SetEnabled(1);
return;
}
# If using reference rasters, make sure the correct number have been selected
if ( toggleDEM.GetValue() ) {
local class GUI_CTRL_LISTBOX list = mainDLG.GetCtrlByID("id_REF_DEM_LIST");
local class GUI_CTRL_EDIT_NUMBER number = mainDLG.GetCtrlByID("id_NUM_REF_DEM");
if ( list.GetCount() != number.GetValue() ) {
PopupMessage("Please select the correct number of DEMs you wish to use as reference");
buttonRun.SetEnabled(1);
return;
}
}
errXML = progressDOC.Parse(progressxmlCODE$);
if ( errXML < 0 ) {
PopupError( errXML ); # pop up an error dialog
mainDLG.Close(1);
Exit( );
}
progressNODE = progressDOC.GetElementByID( "VERIFY_DIALOG" );
if ( progressNODE == 0 ) {
PopupMessage( "Could not find progress dialog node in XML document" );
mainDLG.Close(1);
Exit();
}
progressDLG.SetXMLNode(progressNODE);
progressDLG.DoModal( mainDLG.GetWidget() );
if (numValidNulls == 0) {
print("There are no recognized null cells in the raster.\n");
PopupMessage("There are no recognized null cells in the raster.");
buttonRun.SetEnabled(1);
return;
}
progressNODE = progressDOC.GetElementByID( "PROGRESS_DIALOG" );
if ( progressNODE == 0 ) {
PopupMessage( "Could not find progress dialog node in XML document" );
mainDLG.Close(1);
Exit();
}
progressDLG.SetXMLNode(progressNODE);
# Set up a timer to record completion time
local class TIMER time;
time.Running = 0;
time.Value = 0;
time.Running = 1;
progressDLG.DoModal( mainDLG.GetWidget() );
time.Running = 0;
printf("Time to complete: %2d minutes %0.2d seconds", floor(time.Value / 60), floor(time.Value % 60));
mainDLG.Close(1);
}
######################################################################################################
# Called when "Canel" is pressed. Closes dialog and exits script.
proc OnCancelPressed( ) {
mainDLG.Close(1);
}
######################################################################################################
# Called when "Input SRTM raster" is pressed. Opens dialog to select raster.
proc OnGetInputSRTMPressed( ) {
local class RASTER R;
GetInputRaster(R);
InputRasterFilename$ = R.$Info.Filename;
InputRasterObjectName$ = R.$Info.Name;
rows = NumLins(R);
cols = NumCols(R);
initType = RastType(R);
# Update the dialog
local class GUI_CTRL_EDIT_STRING text = mainDLG.GetCtrlByID("id_INPUT_SRTM_RASTER");
text.SetValue(InputRasterFilename$ + " / " + InputRasterObjectName$, 0);
mainDLG.GetCtrlByID("id_GET_OUTPUT_SRTM_BUTTON").SetEnabled(1);
CloseRaster(R);
}
######################################################################################################
# Called when "Output SRTM raster" is pressed. Opens dialog to select/create output raster.
proc OnGetOutputSRTMPressed( ) {
local class RASTER Rin, Rout;
# Get the raster to save result too
GetOutputRaster(Rout, rows, cols, initType);
OutputRasterFilename$ = Rout.$Info.Filename;
OutputRasterObjectName$ = Rout.$Info.Name;
CloseRaster(Rout);
# Update the dialog
local class GUI_CTRL_EDIT_STRING text = mainDLG.GetCtrlByID("id_OUTPUT_SRTM_RASTER");
text.SetValue(OutputRasterFilename$ + " / " + OutputRasterObjectName$, 0);
mainDLG.GetCtrlByID("id_RUN_BUTTON").SetEnabled(1);
}
######################################################################################################
# Called when toggle is changed. Updates dialog settings.
proc OnUseSingleCellToggleChanged( ) {
local class GUI_CTRL_TOGGLEBUTTON toggle = mainDLG.GetCtrlByID("id_FILL_SINGLE_NULLS_TOGGLE");
mainDLG.GetCtrlByID("id_NUM_REQUIRED_SURROUNDING_LABEL").SetEnabled(toggle.GetValue());
mainDLG.GetCtrlByID("id_NUM_REQUIRED_SURROUNDING").SetEnabled(toggle.GetValue());
}
######################################################################################################
# Called when toggle is changed. Updates dialog settings.
proc OnUseRefDemToggleChanged( ) {
local class GUI_CTRL_TOGGLEBUTTON toggle = mainDLG.GetCtrlByID("id_FILL_VOIDS_WITH_REF_DEM_TOGGLE");
local class GUI_CTRL_EDIT_NUMBER number = mainDLG.GetCtrlByID("id_NUM_REF_DEM");
local class GUI_CTRL_LISTBOX list = mainDLG.GetCtrlByID("id_REF_DEM_LIST");
# Manage the listbox
mainDLG.GetCtrlByID("id_NUM_REF_DEM_LABEL").SetEnabled(toggle.GetValue());
number.SetEnabled(toggle.GetValue());
number.ClearValue(0);
list.DeleteAllItems();
mainDLG.GetCtrlByID("id_GET_REF_DEM_BUTTON").SetEnabled(0);
# Manage the sub-procedure settings for using reference rasters
if (toggle.GetValue() == 0) {
local class GUI_CTRL_TOGGLEBUTTON threshToggle = mainDLG.GetCtrlByID("id_USE_THRESHOLD_TOGGLE");
threshToggle.SetValue(0,1);
threshToggle.SetEnabled(0);
local class GUI_CTRL_TOGGLEBUTTON bufferToggle = mainDLG.GetCtrlByID("id_USE_BUFFER_ZONE_TOGGLE");
bufferToggle.SetValue(0,1);
bufferToggle.SetEnabled(0);
}
}
######################################################################################################
# Called when "Get DEMs..." is pressed. Allows user to select the number of reference DEMs that they
# specified in the number box.
proc OnGetRefDEMRastersPressed( ) {
local class GUI_CTRL_EDIT_NUMBER number = mainDLG.GetCtrlByID("id_NUM_REF_DEM");
local class GUI_CTRL_LISTBOX list = mainDLG.GetCtrlByID("id_REF_DEM_LIST");
list.DeleteAllItems();
RefDEMFilenames$.Clear();
RefDEMObjectNames$.Clear();
local class RASTER R;
local numeric i;
local numeric numRefDems = number.GetValue();
# For the number of reference DEMs the user specified they want to use
for i = 0 to numRefDems-1 {
GetInputRaster(R);
RefDEMFilenames$.AddToEnd(R.$Info.Filename);
RefDEMObjectNames$.AddToEnd(R.$Info.Name);
CloseRaster(R);
list.AddItem(RefDEMFilenames$[i] + " / " + RefDEMObjectNames$[i]);
}
# Update sub-procedures
local class GUI_CTRL_TOGGLEBUTTON threshToggle = mainDLG.GetCtrlByID("id_USE_THRESHOLD_TOGGLE");
threshToggle.SetEnabled(1);
threshToggle.SetValue(0,1);
local class GUI_CTRL_TOGGLEBUTTON bufferToggle = mainDLG.GetCtrlByID("id_USE_BUFFER_ZONE_TOGGLE");
bufferToggle.SetEnabled(1);
bufferToggle.SetValue(0,1);
}
######################################################################################################
# Called after user changes value of the number of reference DEMs they wish to use.
proc OnNumDEMChanged ( ) {
local class GUI_CTRL_EDIT_NUMBER number = mainDLG.GetCtrlByID("id_NUM_REF_DEM");
mainDLG.GetCtrlByID("id_GET_REF_DEM_BUTTON").SetEnabled(1);
local class GUI_CTRL_LISTBOX list = mainDLG.GetCtrlByID("id_REF_DEM_LIST");
list.DeleteAllItems();
RefDEMFilenames$.Clear();
RefDEMObjectNames$.Clear();
# Call procedure to start selecting reference DEMs
OnGetRefDEMRastersPressed();
}
######################################################################################################
# Called when toggle is changed. Updates dialog settings.
proc OnUseThresholdToggleChanged( ) {
local class GUI_CTRL_TOGGLEBUTTON toggle = mainDLG.GetCtrlByID("id_USE_THRESHOLD_TOGGLE");
mainDLG.GetCtrlByID("id_THRESHOLD_VALUE_LABEL").SetEnabled(toggle.GetValue());
mainDLG.GetCtrlByID("id_THRESHOLD_VALUE").SetEnabled(toggle.GetValue());
}
######################################################################################################
# Called when toggle is changed. Updates dialog settings.
proc OnUseBufferZoneToggleChanged( ) {
local class GUI_CTRL_TOGGLEBUTTON toggle = mainDLG.GetCtrlByID("id_USE_BUFFER_ZONE_TOGGLE");
mainDLG.GetCtrlByID("id_NUM_BUFFER_ZONE_LABEL").SetEnabled(toggle.GetValue());
mainDLG.GetCtrlByID("id_NUM_BUFFER_ZONE").SetEnabled(toggle.GetValue());
}
######################################################################################################
# Called when "Project File..." is pressed. Opens dialog to select project file in which to save the
# temporary working vectors and rasters.
proc OnGetTempFilePressed( ) {
local class GUI_CTRL_EDIT_STRING text = mainDLG.GetCtrlByID("id_TEMP_FILE");
TempSRTMfile = GetOutputFileName(_context.ScriptDir, "Select project file to save temp files to:", "rvc");
# Check if the user pressed "Cancel" in the dialog
if (TempSRTMfile == _context.ScriptDir) {
local class GUI_CTRL_TOGGLEBUTTON toggle = mainDLG.GetCtrlByID("id_SAVE_TEMP_FILE_TOGGLE");
mainDLG.GetCtrlByID("id_SAVE_TEMP_FILE_BUTTON").SetEnabled(0);
toggle.SetValue(0, 0);
text.ClearValue(0);
text.SetEnabled(0);
return;
}
# Delete the file if it exists. The user would have recieved a "Overwrite this file?" warning.
if (fexists(TempSRTMfile)) {
DeleteFile(TempSRTMfile);
}
CreateProjectFile(TempSRTMfile, "Created by SRTM Void Filling Script to save temporary files");
text.SetValue(TempSRTMfile,0);
}
######################################################################################################
# Called when toggle is changed. Updates dialog settings.
proc OnSaveTempFileToggleChanged( ) {
local class GUI_CTRL_TOGGLEBUTTON toggle = mainDLG.GetCtrlByID("id_SAVE_TEMP_FILE_TOGGLE");
if (toggle.GetValue()) {
OnGetTempFilePressed();
}
mainDLG.GetCtrlByID("id_SAVE_TEMP_FILE_BUTTON").SetEnabled(toggle.GetValue());
mainDLG.GetCtrlByID("id_TEMP_FILE").SetEnabled(toggle.GetValue());
}
######################################################################################################
# Procedure called to set up the dialog using the XML code specified in this script. Dialog is made
# using a Modal function call. When the dialog closes, the script will end.
proc OnCreateDialog( ) {
errXML = mainDOC.Parse(mainxmlCODE$);
if ( errXML < 0 ) {
PopupError( errXML ); # pop up an error dialog
Exit( );
}
mainNODE = mainDOC.GetElementByID( "MAIN_DIALOG" );
if ( mainNODE == 0 ) {
PopupMessage( "Could not find main dialog node in XML document" );
Exit();
}
mainDLG.SetXMLNode(mainNODE);
errDLG = mainDLG.DoModal();
if ( errDLG == -1 ) {
Exit();
}
}
######################################################################################################