######################################################################################################
#
# SampleCatchments.sml (Stand Alone)
#
######################################################################################################
#
# Author:
# Dave Breitwisch, MicroImages, Inc.
#
# Started:
# 19 September 2005
#
# Updated:
# 10 April 2008
#
# Requires: TNTmips2005:72 or later
#
numeric version = round(_context.Version * 10);
if (version < 72) {
PopupMessage("This version requires that you are using TNTmips2006:72 or later");
Exit();
}
#
######################################################################################################
#
# Documentation: SampleCatchments.sml
#
# Purpose:
#
# The purpose of this Stand Alone SML script is to process multiple vector points representing stream or
# stream sediment samples to compute the upstream catchment (basin) area that drains to each point. Script inputs
# include an elevation raster (DEM) of the area and a vector object containing the sample points with attached attributes.
# The script uses data structures and methods associated with the TNTmips Watershed process.
#
# Script outputs include:
# 1) a vector object containing standard flowpaths determined from the elevation model
# 2) a vector object containing the sample points relocated to the lie on the nearest flowpath (if points
# were within the user-specified tolerance distance)
# 3) a vector object containing one catchment area for each relocated sample point. Each catchment is
# represented by one or more polygon elements. Sample point attributes and computed catchment attributes
# are attached to all associated catchment polygons. Multiple sample points within a single watershed
# result in multiple catchment areas. Catchment attributes track the identity of the immediately adjacent
# upstream samples/catchments and the total area of all contributing catchments.
# 4) optional vector objects containing
# a) upstream flowpaths from each sample point
# b) downstream flowpaths from each sample point
# c) standard catchment basins for the entire area
# d) standard ridges for the entire area
#
# The results from this process can then be analyzed based on the geochemical or other sample data attached
# to the sample points.
#
# The DEM used as input can be either a standard DEM or an adjusted (depressionless) DEM produced by processing
# a standard DEM using the TNTmips Watershed process. For most efficient processing of many sample points,
# first use the Watershed process to create a depressionless DEM and to determine flowpath and basin parameters
# that create a standard flowpath network that covers the input sample set with the minimum density and complexity.
# If using a depressionless DEM as input, turn off the Fill All Depressions toggle on the main dialog window.
#
# Main Dialog:
#
# Input DEM... - Press this button to open up a selection dialog and select the DEM Raster to be used for the watershed process.
#
# Input Vector... - Press this button to open a selection dialog and select the Vector that contains a point database with
# the geochemical samples. The table should be named "TABLE" and the unique sample ID field should be labeled "SAMP_NO".
# Should the database have different names, the script can be altered on line 72 and 73 to match the data if necessary.
#
# Sample Table Name - When selecting an input vector, it will prompt for the user to select the database table for that vector
# that contains the stream sediment samples
#
# Unique ID Field - Prompted for as above for user to specify the field name that contains the unique sample ID values.
#
# Watershed Options -
#
# General - This page allows the user to specify which object will be generated by the script. Those that are grayed out, are
# are either required or not available at the time this script is published.
#
# Mask - This option lets the user select a mask raster for watershed process to use.
#
# Flow Path and Basin - These settings are the same as those in the Watershed process in TNTmips. These allow the user
# to specify how the standard watershed objects are created. For more details on how to use these settings, see the
# "Modeling Watershed Geomorphology" tutorial booklet that is available with your TNTmips installation materials or on the
# MicroImages website.
#
# Snap to closest flowpath toggle - This option is currently always active. Although the user can still specify the maximum
# distance in which to search for the closest standard flowpath line. If there is not a standard flowpath line within the
# specified distance, the user is given a warning and the point is not added to the resulting vector and analysis results.
#
# Snap to downstream node toggle - This option is used to move sample points to the downstream node of the closest line if
# it's within the specified distance. The purpose of this option is to correct a situation in which a sample point may be
# taken downstream of the intersection of two streams in the field, but the "flow path" generated from the DEM Raster ends
# up with a model where the point will be moved to one of the upstream lines above the intersection. As a result of this
# situation, only one upstream line is sampled. By using this option, the script will search within the specified search
# distance downstream for an intersection of another line. If found, the sample point moved to this node and will be marked
# in the database as snapped to a particular node ("SnappedToNode" field). The watershed process and the other computations
# of this script will then sample everything upstream of this node.
#
# Output Project File... - Press this button to open a selection dialog and create a new project file for the output objects
# of the script to be created in.
#
# Output Summary File... - Press this button to open a selection dialog and select a new text file. A summary of certain
# steps, warnings and corrections are printed to both the SML console window and this file. This is an extra tool for the
# user to look through to help evaluate the results of the script.
#
#
# Output RVC Objects Created By The Script:
#
# The following objects (alphabetical order) are created in the project file created when the "Output Project File..." button was pressed.
#
# Result_Vector - This vector is created with a copy of the point database table that is in the Vector selected when the
# "Input Vector..." button was pressed. After the STDFLOWPATH Vector has been created, each sample point is moved to the
# closest line and then to the center of the DEM Raster cell on that line (all flow paths are computed to go through the
# center of the cells). This point is added to the vector and attached to the appropriate record in the table. If a point
# is to be added in the same location (center of same DEM Raster cell), the sample ID is added to a record (see description
# of "DupSample" and "NumDuplicates" fields below) of the point that already exists and the point is not added to the vector.
#
# Additional fields are added to the copied Vector point database table (ResultSampleTBL):
#
# "TABLE" - This is the current name of the table that is to be copied from the "Input Vector", and the name of the table
# once it is copied to the Result_Vector.
#
# "ModEasting" - This is the Easting map coordinate of the sample point after it has been moved to the center of the DEM
# Raster cell of the nearest flow path.
#
# "ModNorthing" - This is the Northing map coordinate of the sample point after it has been moved to the center of the DEM
# Raster cell of the nearest flow path.
#
# "Displaced" - This is the distance in meters in which the point was moved from it's original location to the new location.
#
# "SnappedToNode" - This is the node element number that the point was moved to if it was in the maximum search distance
# specified by the user. This value will only be set if the "Snap to downstream node" toggle has been enabled. The value
# of (-1) is the default value for invalid.
#
# "Horton" - The Horton stream order value of the flowpath line that the sample point was moved to.
#
# "Shreve" - The Shreve stream order value of the flowpath line that the sample point was moved to.
#
# "Strahler" - The Strahler stream order value of the flowpath line that the sample point was moved to.
#
# "Scheidegger" - The Scheidegger stream order value of the flowpath line that the sample point was moved to.
#
# "TotalUpSamples" - This is the total number of samples upstream of the current sample point based on the flow path created.
#
# "TotalUpChanLen" - This is the sum of all flow path line segments that are considered upstream of the sample point. This
# also includes the length from the sample point to the nearest upstream node. This value is recorded in kilometers (km).
#
# "TotalArea_sq_km" - This is the area of the total upstream catch basin for the sample ID attached to the record.
#
# The following new table (ResultPointToPointTBL) and fields are created in the Vector point database:
#
# "PointToPoint" - This is the current name of the point database table that is created to show relationships between the
# sample points in the vector.
#
# "OnLine" - This is the line element number in the STDFLOWPATH Vector that the sample point has been moved to. If a sample
# point was moved to the location of a node, the element number of the line downstream of that node is used.
#
# "NumDuplicates" - This specifies the number of sample points that were to be moved to the same location. Because the
# duplicate points were not added to the vector, they are listed in the "DupSample#" fields.
#
# "DupSample#" - This specifies the sample ID of the duplicate point. The '#' in the field name represents incrementing
# numbers (1 to the maximum number of duplicates found for a single point), so each duplicate sample ID for a particular
# point has its own field. These fields are added by the script as necessary.
#
# "DownSample" - This is the sample ID of the next sample point downstream. A value of (-1) represents no samples
# downstream based on the flow path created.
#
# "TotalUpSamples" - This is the total number of samples upstream of the current sample point based on the flow path created.
#
# "NumUpSamples" - This is the number of samples that are upstream of the current sample point based on the flowpath
# created. This only identifies the number that are immediately reachable by tracing up the flow path without tracing
# through an additional sample point.
#
# "UpSample#" - This is the sample ID of an upstream sample point that is immediately reachable by tracing up the flow path
# without tracing through any additional sample points. The '#' in the field name represents incrementing numbers (1 to the
# maximum number of upstream sample points for any point), so each upstream sample point has it's own field for that point.
#
# STDBASIN (Optional) - This is a vector object that is created by the standard watershed process. The basins computed in this vector
# are based on seed points generated by the standard process. This is copied to the output project file by the script, but
# otherwise not used in the additional script computations.
#
# STDFLOWPATH - This is a vector object that is created by the standard watershed process. The flow paths computed in this
# vector are based on seed points generated by the standard process.
#
# The following new table (STDFlowPathTBL) and fields are created in the Vector line database:
#
# "LineToLine" - This is the current name of the Vector line database table that is created to show up and down stream
# relationships between the line segments in the vector.
#
# "DownNode" - This is the node element number that is considered the downstream endpoint of the line segment.
#
# "DownLine" - This is the element number of the line that is attached at the "DownNode" and is considered downstream of the
# current line segment.
#
# "UpNode" - This is the element number of the node that is considered the upstream endpoint of the line segment.
#
# "NumUpLines" - This number represents the number of lines attached to the "UpNode" of the line segment. These lines are
# all considered to be upstream of the current line. The current line is excluded in this number of lines attached to the node.
#
# "UpLine#" - This is the element number of the line that is attached to the "UpNode" and is considered upstream of the
# current line. The '#' in the field name represents incrementing numbers from 1 to 7 (the maximum number of possible
# upstream lines attached at the same node), so each upstream line number has it's own field.
#
# "DownSample" - This is the sample ID of the sample point that is located furthest downstream on the line.
#
# "UpSample" - This is the sample ID of the sample point that is located furthest upstream on the line.
#
# "UpChanLen" - This is the sum of lengths of all line segments that are considered upstream of the current line. This
# number is computed in kilometers (km). The length of the current line is excluded from this number.
#
# STDRIDGE (Optional) - This is a vector object that is created by the standard watershed process. The ridge lines computed in this
# vector are based on the standard process. This is copied to the output project file by the script but otherwise not used
# in the additional script computations.
#
# USDBASIN - This is a vector object that is created by the watershed process. The basins are generated using the points
# that were modified and added to the Result_Vector as seed points for the watershed process.
#
# Additional fields are added to the copied Vector polygon sediment database table (UBasinSampleTBL):
#
# "TotalArea_sq_km" - This is the area of the total upstream catch basin for the sample ID attached to the record. This
# includes all upstream basins and the basin for the sample point related by the sample ID.
#
# The following new table (UBasinPolyToPolyTBL) and fields are created in the Vector polygon database:
#
# "PolyToPoly" - This is the current name of the polygon database table that is created to show relationships between the
# user defined watershed basins (polygons) in the vector.
#
# "SAMP_NO" - This is the sample ID of the sample point that was used to create the polygon. This field name should remain
# The same as in the original database table that was copied to the Result_Vector.
#
# "NumPolyPerSamp" - This is the number of Vector polygons that represent the basin of a single sample ID. These polygons
# are all attached to this record.
#
# "SampleBasinArea" - This is the sum of the areas for all polygons that represent the basin of a single sample ID. This
# does not include basins that are created upstream due to additional sample points.
#
# "TotalSampIncl" - This is the number of sample points that are included in the total upstream catch basin for the sample
# point. This includes the sample point attached to this record.
#
# "TotalArea_sq_km" - This is the area of the total upstream catch basin for the sample ID attached to the record. This
# includes all upstream basins and the basin for the sample point related by the sample ID.
#
# "TotalUpBasins" - This is the number of basins that make up the total catch basin for the sample point related by the sample ID.
#
# "NumUpBasins" - This is the number of basins that are upstream basins that are adjacent to the current basin. This is
# the same value as "NumUpSamples" in the "PointToPoint" table in the Result_Vector.
#
#
# USDFLOWPATH (Optional) - This is a vector object that is created by the watershed process. The flow paths are generated using the
# points that were modified and added to the Result_Vector as seed points for the watershed process. This is copied to the
# output project file by the script but otherwise not used in the additional script computations.
#
#
# USD_UPFLOWPATH (Optional) - This is the STDFLOWPATH vector that has been clipped using the USDBASIN vector. Therefore it is considered
# an upstream flowpath vector of all the user defined sample points.
#
######################################################################################################
######################################################################################################
# Global Variable Declarations
######################################################################################################
numeric DEBUG = 1;
class GUI_DLG mainDLG;
class STATUSDIALOG statusD;
class STATUSCONTEXT statusC;
numeric errXML, errDLG;
string xmlCODE$;
string InputDEMfilename = "", InputVECTORfilename = "", InputMaskfilename;
string InputDEMobjectName, InputVECTORobjectName, InputMaskobjectName;
string OutputProjectFile = _context.ScriptDir;
string OutputSummaryFile = _context.ScriptDir;
class FILE summaryFile;
string message;
class WATERSHED watershed;
class RASTER SampleDEM;
class VECTOR SampleVECT;
class VECTOR WarpVECT;
class VECTOR ResultVECT;
class VECTOR UBasinVECT;
class TRANSPARM warptransparm;
##########################################################
# Define the table's names, descriptions, and fields here:
##########################################################
#
# Names must be kept to 15 charactors or less. Numbers
# added to the field name will count against this limit.
#
##########################################################
######################
# Result_Vector Tables
######################
# ResultSampleTBL = "TABLE"
string sampleTblName = "TABLE";
string sampleIDFldName = "SAMP_NO"; # Sample point unique ID
string modEastFldName = "ModEasting"; # Modified Easting value of the moved point
string modNorthFldName = "ModNorthing"; # Modified Northing value of the moved point
string displacedFldName = "Displaced"; # Distance the point was moved from its original location
string hortonFldName = "Horton"; # "Horton" Stream order
string shreveFldName = "Shreve"; # "Shreve" Stream order
string strahlerFldName = "Strahler"; # "Stahler" Stream order
string scheideggerFldName = "Scheidegger"; # "Scheidegger" Stream order
string snapToDownNodeFldName = "SnappedToNode"; # Node the sample was snapped to
string totalUpChanLenFldName = "TotalUpChanLen"; # Length of the upstream channel
# ResultPointToPointTBL => Table created to show the relationships between the sample points
string PointToPointTblName = "PointToPoint";
string PointToPointTblDesc = "Point to Point Relationships";
string onLineFldName = "OnLine"; # Line number that the sample point is on
string numDuplFldName = "NumDuplicates"; # Number of duplicate samples in the same DEM cell
string dupSampleIDFldName = "DupSample";# + Number # Duplicate sample point ID, added by script if necessary
string downSampleFldName = "DownSample"; # First downstream sample encountered
string totalUpSamplesFldName = "TotalUpSamples"; # Total number of samples upstream
string numUpSamplesFldName = "NumUpSamples"; # Number of immediate upstream samples
string upSampleFldName = "UpSample";# + Number # Immediate upstream sample
######################
# STDFLOWPATH Table
######################
# STDFlowPathTBL => Table created to show the relationships between the lines in the flowpath vector that is created
string LineToLineTblName = "LineToLine";
string LineToLineTblDesc = "Line to Line Flowpath Relationships";
string downNodeFldName = "DownNode"; # The end node of the line that is downstream
string downLineFldName = "DownLine"; # The connected line that is downstream
string upNodeFldName = "UpNode"; # The end node of the line that is upstream
string numUpLinesFldName = "NumUpLines"; # The number of connected lines that are uptream of the UpNode
string upLineFldName = "UpLine";# + Number # A line connected to the upnode that is upstream
string mostDownSampleFldName = "DownSample"; # The most downstream sample on the line segment
string mostUpSampleFldName = "UpSample"; # The most upstream sample on the line segment
string totalUpChanLengthFldName = "UpChanLen"; # The total length of all lines upstream of the line segment (self excluded)
######################
#USDBASIN Table
######################
# UBasinPolyToPoly => Table created to show the relationships between the user-defined basins
string PolyToPolyTblName = "PolyToPoly";
string PolyToPolyTblDesc = "Polygon to Polygon Basin Relationships";
#string sampleIDFldName = "SAMP_NO"; # Sample ID realated to the basin polygon. SAME AS ResultSampleTBL
string numPolysPerSampleFldName = "NumPolyPerSamp"; # Number of vector polygons related to that sample
string areaSamplePolysFldName = "SampBasinArea"; # Total of all vector polygons related to that sample
string totalSamplesIncludedFldName = "TotalSampIncl"; # Number of Samples included (self included)
string totalSampleBasinAreaFldName = "TotalArea_sq_km"; # Total area of all polygons upstream of sample point (self included)
string totalSampleBasinsFldName = "TotalUpBasins"; # Total number of sample basins for sample point (self included)
string numUpBasinsFldName = "NumUpBasins"; # Number of basins upstream (self excluded)
######################################################################################################
######################################################################################################
# Dialog XML Code
######################################################################################################
xmlCODE$ = '
';
######################################################################################################
######################################################################################################
######################################################################################################
# Start of the Script
######################################################################################################
clear();
proc OnCreateDialog(); #Function Prototype
OnCreateDialog();
######################################################################################################
######################################################################################################
func string SampleIDThatIsNotUnique ( ) {
##################################################################################################
# Checking to make sure that all the sample IDs in the point vector table are unique
statusC.Message = "Verifying all sample IDs are unique";
##################################################################################################
local class VECTOR UniqueCheckVector;
OpenVector(UniqueCheckVector, InputVECTORfilename, InputVECTORobjectName);
local class DATABASE pntDBase = OpenVectorPointDatabase(UniqueCheckVector);
local class DBTABLEINFO tblInfo = DatabaseGetTableInfo(pntDBase, sampleTblName);
local numeric numTotalRecords = tblInfo.NumRecords;
numeric uniquehash[];
local numeric index;
local string key;
for index = 1 to numTotalRecords {
key = TableReadFieldStr(tblInfo, sampleIDFldName, index);
if ( uniquehash.Exists(key) ) {
uniquehash.Clear();
print(i, key, sampleIDFldName);
return key;
} else {
uniquehash[key] = 1;
}
}
uniquehash.Clear();
return "";
}
######################################################################################################
######################################################################################################
proc StartWatershedProcess ( ) {
##################################################################################################
# Initialize the watershed objects and create the output objects for the script
statusC.Message = "Step 1 of 18: Computing standard Watershed objects";
##################################################################################################
watershed = WatershedInit(InputDEMfilename, InputDEMobjectName);
WatershedSetBasin(watershed, mainDLG.GetCtrlByID("id_BASIN").GetValueNum());
WatershedSetInlet(watershed, mainDLG.GetCtrlByID("id_INLET").GetValueNum());
WatershedSetBranch(watershed, mainDLG.GetCtrlByID("id_BRANCH").GetValueNum());
WatershedSetOutlet(watershed, mainDLG.GetCtrlByID("id_OUTLET").GetValueNum());
if (mainDLG.GetCtrlByID("id_SEP_VALLEY_POLYS").GetValueNum()) {
WatershedSetValleySeparation(watershed, mainDLG.GetCtrlByID("id_BASIN").GetValueNum());
}
if (mainDLG.GetCtrlByID("id_MASK_TOGGLE").GetValueNum()) {
if (ObjectExists(InputMaskfilename, InputMaskobjectName, "Raster")) {
WatershedSetMask(watershed, InputMaskfilename, InputMaskobjectName);
} else {
message = "##########################################################################################\n" +
"Waning: The the Mask Raster selected for the watershed process is invalid.\n" +
" A full watershed analysis will be done instead.\n" +
"##########################################################################################\n\n";
print(message);
fwritestring(summaryFile, message);
}
}
# Set up the watershed object creation parameters
local string watershedparms$ = "";
if ( mainDLG.GetCtrlByID("id_FILL_ALL_DEP").GetValueNum()) {
watershedparms$ += "FillAllDepressions,";
}
if ( mainDLG.GetCtrlByID("id_COMP_STD_BASINS").GetValueNum()) {
watershedparms$ += "Basin,";
}
if ( mainDLG.GetCtrlByID("id_COMP_STD_RIDGES").GetValueNum()) {
watershedparms$ += "Ridge,";
}
# if ( mainDLG.GetCtrlByID("id_COMP_STD_DBASE").GetValueNum()) {
# watershedparms$ += "Database,";
# }
watershedparms$ += "FlowPath";
print(watershedparms$);
WatershedCompute(watershed, watershedparms$);
#array numeric xvals[0];
#array numeric yvals[0];
#WatershedComputeElements(watershed, xvals, yvals, 0, "FlowPath");
# Copy Standard Watershed Objects
local string ObjectFilename, ObjectName;
local numeric ObjNumber;
WatershedGetObject(watershed, "VectorFlowPath", ObjectFilename, ObjectName);
ObjNumber = ObjectNumber(ObjectFilename, ObjectName, "VECTOR");
CopyObject(ObjectFilename, ObjNumber, OutputProjectFile);
if ( mainDLG.GetCtrlByID("id_COMP_STD_BASINS").GetValueNum()) {
WatershedGetObject(watershed, "VectorBasin", ObjectFilename, ObjectName);
ObjNumber = ObjectNumber(ObjectFilename, ObjectName, "VECTOR");
CopyObject(ObjectFilename, ObjNumber, OutputProjectFile);
}
if ( mainDLG.GetCtrlByID("id_COMP_STD_RIDGES").GetValueNum()) {
WatershedGetObject(watershed, "VectorRidge", ObjectFilename, ObjectName);
ObjNumber = ObjectNumber(ObjectFilename, ObjectName, "VECTOR");
CopyObject(ObjectFilename, ObjNumber, OutputProjectFile);
}
# Create new vector to contain the result of the script
OpenRaster(SampleDEM, InputDEMfilename, InputDEMobjectName);
OpenVector(SampleVECT, InputVECTORfilename, InputVECTORobjectName);
local class GEOREF DEMgeoref = GetLastUsedGeorefObject(SampleDEM);
local class GEOREF VECTgeoref = GetLastUsedGeorefObject(SampleVECT);
warptransparm.InputCoordRefSys = VECTgeoref.CoordRefSys;
warptransparm.OutputCoordRefSys = DEMgeoref.CoordRefSys;
CreateVector(WarpVECT, OutputProjectFile, "WarpVect", "Warped Input Vector"); #XXX
VectorWarp(SampleVECT, WarpVECT, warptransparm, 10);
# Compute the extents of the WarpVECTOR to use in initializing the results vector.
# VectorValidate(WarpVECT);
local class REGION2D extents = GetObjectExtentsRegion(WarpVECT, GetLastUsedGeorefObject(WarpVECT));
CreateVector(ResultVECT, OutputProjectFile, "Result Vector", "Resulting Vector created by SML Watershed Script", "VectorToolkit", extents.Extents);
# VectorCopyElements(WarpVECT, ResultVECT, "RemoveTables", 0, 0, 0);
CopySubobjects(WarpVECT, ResultVECT, "GEOREF");
local class GEOREF WARPgeoref = GetLastUsedGeorefObject(WarpVECT);
# Open the standard flowpath vector
local class VECTOR STDFlowPath;
OpenVector(STDFlowPath, OutputProjectFile, "STDFLOWPATH");
##################################################################################################
# Setup the Database tables used during the script
statusC.Message = "Step 2 of 18: Creating new database tables";
##################################################################################################
# Copy the sample database table from original vector
local class DATABASE SampleDB = OpenVectorPointDatabase(WarpVECT);
local class DBTABLEINFO SampleTBL = DatabaseGetTableInfo(SampleDB, sampleTblName);
local class DATABASE ResultDB = OpenVectorPointDatabase(ResultVECT);
if (TableCopy(SampleDB, SampleTBL, ResultDB) < 0) {
message = "Error occured while trying to copy the sample table.\n" + "Exiting the script now... \n";
print(message);
fwritestring(summaryFile, message);
CloseRaster(SampleDEM);
CloseVector(SampleVECT);
CloseVector(WarpVECT);
CloseVector(ResultVECT);
Exit();
}
# Add additional fields to the sample table to note movement in the sample locations
local class DBTABLEINFO ResultSampleTBL = DatabaseGetTableInfo(ResultDB, sampleTblName);
TableAddFieldFloat(ResultSampleTBL, modEastFldName, 15, 4); # Modified Easting position of the sample point
TableAddFieldFloat(ResultSampleTBL, modNorthFldName, 15, 4); # Modified Northing position of the sample point
TableAddFieldFloat(ResultSampleTBL, displacedFldName, 10, 3); # Displacement in meters from the original to modified position
TableAddFieldInteger(ResultSampleTBL, snapToDownNodeFldName, 10); # Specifies the downstream node that it was snapped to if the node
# was within the distance and on the closests line
TableAddFieldInteger(ResultSampleTBL, hortonFldName, 7); # "Horton" stream order of the sample
TableAddFieldInteger(ResultSampleTBL, shreveFldName, 7); # "Shreve" stream order of the sample
TableAddFieldInteger(ResultSampleTBL, strahlerFldName, 9); # "Strahler" stream order of the sample
TableAddFieldInteger(ResultSampleTBL, scheideggerFldName, 12); # "Scheidegger" stream order of the sample
TableAddFieldInteger(ResultSampleTBL, totalUpSamplesFldName, 10); # Total number of upstream sample points
TableAddFieldFloat(ResultSampleTBL, totalUpChanLenFldName, 10, 3); # The total upstream channel length being sampled in kilometers
TableAddFieldFloat(ResultSampleTBL, totalSampleBasinAreaFldName, 15, 4); # The total area of the upstream catchment basin
# Create a table to show relationship between the sample points
local class DBTABLEINFO ResultPointToPointTBL;
ResultPointToPointTBL = TableCreate(ResultDB, PointToPointTblName, PointToPointTblDesc);
TableAddFieldInteger(ResultPointToPointTBL, onLineFldName, 10); # Line number that the sample point is on
TableAddFieldInteger(ResultPointToPointTBL, numDuplFldName, 13); # Number of duplicate samples in the same DEM cell
#TableAddFieldString(ResultPointToPointTBL, dupSampleIDFldName + #, 12); # Duplicate sample point ID, added by script if necessary
TableAddFieldString(ResultPointToPointTBL, downSampleFldName, 11); # First downstream sample encountered
TableAddFieldInteger(ResultPointToPointTBL, totalUpSamplesFldName, 10); # Total number of samples upstream
TableAddFieldInteger(ResultPointToPointTBL, numUpSamplesFldName, 10); # Number of immediate upstream samples
TableAddFieldString(ResultPointToPointTBL, upSampleFldName + "1", 11); # Immediate upstream sample
TableAddFieldString(ResultPointToPointTBL, upSampleFldName + "2", 11); # Additional immediate upstream sample
#TableAddFieldString(ResultPointToPointTBL, upSampleFldName + "#", 11); # Additional added by script if necessary
# Create a table to show the relationships between the flowpath lines created by the watershed
# process to help determine the direction of the flowpath. Also shows the most upstream and
# downstream sample points on each line.
local class DATABASE STDFlowPathDB = OpenVectorLineDatabase(STDFlowPath);
local class DBTABLEINFO STDFlowPathTBL;
STDFlowPathTBL = TableCreate(STDFlowPathDB, LineToLineTblName, LineToLineTblDesc);
TableAddFieldInteger(STDFlowPathTBL, downNodeFldName, 10); # The end node of the line that is downstream
TableAddFieldInteger(STDFlowPathTBL, downLineFldName, 10); # The connected line that is downstream
TableAddFieldInteger(STDFlowPathTBL, upNodeFldName, 10); # The end node of the line that is upstream
TableAddFieldInteger(STDFlowPathTBL, numUpLinesFldName, 10); # The number of connected lines that are uptream of the UpNode
TableAddFieldInteger(STDFlowPathTBL, upLineFldName + "1", 10); # A line connected to the upnode that is upstream
TableAddFieldInteger(STDFlowPathTBL, upLineFldName + "2", 10); # A line connected to the upnode that is upstream
TableAddFieldInteger(STDFlowPathTBL, upLineFldName + "3", 10); # A line connected to the upnode that is upstream
TableAddFieldInteger(STDFlowPathTBL, upLineFldName + "4", 10); # A line connected to the upnode that is upstream
TableAddFieldInteger(STDFlowPathTBL, upLineFldName + "5", 10); # A line connected to the upnode that is upstream
TableAddFieldInteger(STDFlowPathTBL, upLineFldName + "6", 10); # A line connected to the upnode that is upstream
TableAddFieldInteger(STDFlowPathTBL, upLineFldName + "7", 10); # A line connected to the upnode that is upstream
TableAddFieldString(STDFlowPathTBL, mostDownSampleFldName, 11); # The most downstream sample on the line segment
TableAddFieldString(STDFlowPathTBL, mostUpSampleFldName, 11); # The most upstream sample on the line segment
TableAddFieldFloat(STDFlowPathTBL, totalUpChanLengthFldName, 15, 3); # The total length of all lines upstream of the line segment (self excluded)
##################################################################################################
# Common variables used in multiple steps
##################################################################################################
local numeric numOriginalSamplePoints;
local numeric numSamplePoints;
local numeric i, j, k, node1, node2;
local class POLYLINE P;
local class POINT2D pt1, pt2;
local array numeric records[100];
local numeric Temp1MapX, Temp1MapY, Temp1ObjX, Temp1ObjY, Temp2MapX, Temp2MapY;
local numeric numFlowPathLines = NumVectorLines(STDFlowPath);
local string tempFieldName;
local string sampleNum$;
local array numeric elements[10];
local class STRINGLIST UpSample;
local string tempStr$;
local string tempID1, tempID2;
local numeric warnCountDuplicates = 0;
local numeric warnCountOutsideSearch = 0;
local numeric warnCountLostPolygons = 0;
local numeric warnCountManAttachPolygons = 0;
local numeric numComputed;
local array numeric lines1[10], lines2[10];
# Get the georeference objects for transformations
local class GEOREF FLOWgeoref = GetLastUsedGeorefObject(STDFlowPath);
local class GEOREF RSLTgeoref = GetLastUsedGeorefObject(ResultVECT);
##################################################################################################
# Determine the upstream and downstream relationship based on the STDFlowPath Vector
statusC.Message = "Step 3 of 18: Computing up and downstream flowpath relationships";
##################################################################################################
# Variables used
local numeric numAttachedLines1, numAttachedLines2, currentShreveValue;
local array upLines[7];
local numeric downLine, upNode, downNode, lineCount, numUpLines;
# For each line in the STDFlowPath Vector
statusC.BarInit(numFlowPathLines, 0);
for i = 1 to numFlowPathLines {
statusC.BarIncrement(1, 0);
# Initialize the variables to invalid
downNode = -1;
downLine = -1;
upNode = -1;
for j = 1 to 7 {
upLines[j] = -1;
}
currentShreveValue = STDFlowPath.line[i].STREAM_ORDER.Shreve;
# Get the line
P = GetVectorLine(STDFlowPath, i);
# Get the first end point and the lines attached to it
pt1 = P.GetVertex(0);
pt1 = ObjectToMap(STDFlowPath, pt1.x, pt1.y, FLOWgeoref);
node1 = FindClosestNode(STDFlowPath, pt1.x, pt1.y, FLOWgeoref);
numAttachedLines1 = GetVectorNodeLineList(STDFlowPath, lines1, node1);
# Get the second end point and the lines attached to it
pt2 = P.GetVertex(P.GetNumPoints() - 1);
pt2 = ObjectToMap(STDFlowPath, pt2.x, pt2.y, FLOWgeoref);
node2 = FindClosestNode(STDFlowPath, pt2.x, pt2.y, FLOWgeoref);
numAttachedLines2 = GetVectorNodeLineList(STDFlowPath, lines2, node2);
# If the only lines attached to both nodes is the line connecting the two nodes
if ( (numAttachedLines1 == 1) and (numAttachedLines2 == 1) ) {
numUpLines = 0;
# Determine upstream based on Z values of the nodes. Tie sets upstream to node1
if (STDFlowPath.Node[node1].Internal.z >= STDFlowPath.Node[node2].Internal.z) {
downNode = node2;
upNode = node1;
} else {
downNode = node1;
upNode = node2;
}
# If node1 has multiple attached lines and node2 has 1
} else if ( (numAttachedLines1 > 1) and (numAttachedLines2 == 1) ) {
# For each line attached to node1
for j = 1 to numAttachedLines1 {
# Don't check the current line
if (lines1[j] != i) {
# Check the Shreve stream order system. Greatest value will be downstream
if (STDFlowPath.line[lines1[j]].STREAM_ORDER.Shreve > currentShreveValue) {
currentShreveValue = STDFlowPath.line[lines1[j]].STREAM_ORDER.Shreve;
downLine = lines1[j];
downNode = node1;
upNode = node2;
}
}
}
# If node1 is not downstream, node2 must be
if (downNode == -1) {
downNode = node2;
upNode = node1;
}
} else if ( (numAttachedLines1 == 1) and (numAttachedLines2 > 1) ) {
# For each line attached to node2
for j = 1 to numAttachedLines2 {
# Don't check the current line
if (lines2[j] != i) {
# Check the Shreve stream order system. Greatest value will be downstream
if (STDFlowPath.line[lines2[j]].STREAM_ORDER.Shreve > currentShreveValue) {
currentShreveValue = STDFlowPath.line[lines2[j]].STREAM_ORDER.Shreve;
downLine = lines2[j];
downNode = node2;
upNode = node1;
}
}
}
# If node2 is not downstream, node1 must be
if (downNode == -1) {
downNode = node1;
upNode = node2;
}
# Both nodes have multiple attached lines
} else {
# For each line attached to node1
for j = 1 to numAttachedLines1 {
# Don't check the current line
if (lines1[j] != i) {
# Check the Shreve stream order system. Greatest value will be downstream
if (STDFlowPath.line[lines1[j]].STREAM_ORDER.Shreve > currentShreveValue) {
currentShreveValue = STDFlowPath.line[lines1[j]].STREAM_ORDER.Shreve;
downLine = lines1[j];
downNode = node1;
upNode = node2;
}
}
}
# For each line attached to node2
for j = 1 to numAttachedLines2 {
# Don't check the current line
if (lines2[j] != i) {
# Check the Shreve stream order system. Greatest value will be downstream
if (STDFlowPath.line[lines2[j]].STREAM_ORDER.Shreve > currentShreveValue) {
currentShreveValue = STDFlowPath.line[lines2[j]].STREAM_ORDER.Shreve;
downLine = lines2[j];
downNode = node2;
upNode = node1;
}
}
}
} # End of finding nodes and attached lines
# If node1 is the down node, then lines2 are the upstream lines
lineCount = 0;
if (downNode == node1) {
numUpLines = numAttachedLines2 - 1;
for k = 1 to numAttachedLines2 {
if (lines2[k] != i) {
lineCount = lineCount + 1;
upLines[lineCount] = lines2[k];
}
}
# If node2 is the down node, then lines1 are the upstream lines
} else if (downNode == node2) {
numUpLines = numAttachedLines1 - 1;
for k = 1 to numAttachedLines1 {
if (lines1[k] != i) {
lineCount = lineCount + 1;
upLines[lineCount] = lines1[k];
}
}
# Print message if unable to determine.... Shouldn't happen though
} else {
message = "Warning: Unable to determine upstream and downstream configuration for line: " + NumToStr(i) + "\n\n";
print(message);
fwritestring(summaryFile, message);
}
# Write record of the line that is downstream
records[1] = TableWriteRecord(STDFlowPathTBL, 0, downNode, downLine, upNode, numUpLines,
upLines[1], upLines[2], upLines[3], upLines[4], upLines[5], upLines[6], upLines[7]);
TableWriteAttachment(STDFlowPathTBL, i, records, 1, "line");
}
##################################################################################################
# Determine the upstream channel lengths
statusC.Message = "Step 4 of 18: Computing upstream channel lengths for each flowpath line";
##################################################################################################
local numeric upLine, currUpLength, tempUpLength, tempLength, tempCount = 0;
numComputed = 0;
# For each line in the STDFlowPath Vector
statusC.BarInit(numFlowPathLines, 0);
for i = 1 to numFlowPathLines {
statusC.BarIncrement(1, 0);
TableReadAttachment(STDFlowPathTBL, i, records, "line");
upLine = TableReadFieldNum(STDFlowPathTBL, upLineFldName + "1", records[1]);
# If the current line does not have an upstream line, set it's Upstream Channel Length to 0;
if (upLine == -1) {
TableWriteField(STDFlowPathTBL, records[1], totalUpChanLengthFldName, 0);
numComputed = numComputed + 1;
# Otherwise set it to invalid stating it needs to be computed
} else {
TableWriteField(STDFlowPathTBL, records[1], totalUpChanLengthFldName, -1);
}
}
# Until all upstream field that can be computed, are computed
statusC.BarInit(numFlowPathLines, 0);
while (numComputed > 0) {
# Reset the counter and go though all the flowpath lines
tempCount += numComputed;
statusC.BarIncrement(numComputed, 0);
numComputed = 0;
for i = 1 to numFlowPathLines {
TableReadAttachment(STDFlowPathTBL, i, records, "line");
currUpLength = TableReadFieldNum(STDFlowPathTBL, totalUpChanLengthFldName, records[1]);
# If the Upstream length has not been computed
if (currUpLength == -1) {
numUpLines = TableReadFieldNum(STDFlowPathTBL, numUpLinesFldName, records[1]);
tempUpLength = 0;
# For each upstream line
for j = 1 to numUpLines {
# Check that it hasn't failed yet
if (tempUpLength != -1) {
# Get one of the Upstream lines
TableReadAttachment(STDFlowPathTBL, i, records, "line");
upLine = TableReadFieldNum(STDFlowPathTBL, upLineFldName + NumToStr(j), records[1]);
# Get it's upstream channel length value
TableReadAttachment(STDFlowPathTBL, upLine, records, "line");
tempLength = TableReadFieldNum(STDFlowPathTBL, totalUpChanLengthFldName, records[1]);
# If it has been computed, add the value and the length of this upstream line to the temp total
if (tempLength != -1) {
tempUpLength = tempUpLength + tempLength;
tempUpLength = tempUpLength + (STDFlowPath.line[upLine].LINESTATS.Length / 1000);
# If it has not been computed, invalidate the total and move on
} else {
tempUpLength = -1;
}
}
}
# If the Upstream Channel length was successfully computed
if (tempUpLength != -1) {
numComputed = numComputed + 1;
TableReadAttachment(STDFlowPathTBL, i, records, "line");
TableWriteField(STDFlowPathTBL, records[1], totalUpChanLengthFldName, tempUpLength);
}
} # End if currUpLength needs to be computed
} # End for i loop
} # End while
##################################################################################################
# Move the sample points to the closest flowpath line. If close enough, move to the the end node
# of that line if the displacement distance is within the user specified maximum value. Once on
# the line, move to the center of that DEM raster cell.
statusC.Message = "Step 5 of 18: Computing best sample point locations based on user supplied seed points";
##################################################################################################
# Determine the number of relevant vector points that will be used
local numeric numOriginalSamplePoints = NumVectorPoints(WarpVECT);
local numeric currNumDups, currMaxNumDups = 0, numPointsAdded = 0;
array numeric x[numOriginalSamplePoints], y[numOriginalSamplePoints];
# Move each point to the closest flowpath and add to result vector
local numeric snappedToNode, pointExists;
local numeric maxSnapNodeDist = mainDLG.GetCtrlByID("id_SNAP_NODE_MAX_DIST").GetValueNum();
local numeric maxSnapLineDist = mainDLG.GetCtrlByID("id_SNAP_LINE_MAX_DIST").GetValueNum();
local numeric pt1Dist, pt2Dist, tempDist;
local numeric OrigMapX, OrigMapY, OrigObjX, OrigObjY;
local numeric MovedMapX, MovedMapY, MovedObjX, MovedObjY;
local numeric VectObjX, VectObjY;
local numeric err, lineNum, numRecords;
local numeric distance, azimuth, elevation; #azimuth and elev are not used
local class DBFIELDINFO beforeField = FieldGetInfoByName(ResultPointToPointTBL, downSampleFldName);
statusC.BarInit(numOriginalSamplePoints, 0);
for i = 1 to numOriginalSamplePoints {
statusC.BarIncrement(1, 0);
OrigObjX = WarpVECT.Point[i].Internal.x; # In object coordinates
OrigObjY = WarpVECT.Point[i].Internal.y;
ObjectToMap(WarpVECT, OrigObjX, OrigObjY, WARPgeoref, OrigMapX, OrigMapY);
# Move to the closed point on the closest flowpath line in object coordinates
lineNum = FindClosestLine(STDFlowPath, OrigMapX, OrigMapY, FLOWgeoref, maxSnapLineDist);
# Ignore the point if the line that it is attached to does not have a stream order record assosiated with it
if (TableReadAttachment(STDFlowPathTBL, lineNum, records, "line") == 0) {
TableReadAttachment(SampleTBL, i, records, "point");
tempID1 = TableReadFieldStr(SampleTBL, sampleIDFldName, records[1]);
warnCountOutsideSearch = warnCountOutsideSearch + 1;
message = "Warning: Point " + tempID1 + " was to be moved to the standard flowpath line " + NumToStr(lineNum) + ", but the stream order for\n" +
" for this line was unable to be determined. Since the script uses the standard\n" +
" flowpath vector to compute up and downstream relationships, this point is being\n" +
" discarded and not used.\n\n";
print(message);
fwritestring(summaryFile, message);
continue;
}
MapToObject(FLOWgeoref, OrigMapX, OrigMapY, STDFlowPath, OrigObjX, OrigObjY);
ClosestPointOnLine(STDFlowPath, lineNum, OrigObjX, OrigObjY, MovedObjX, MovedObjY);
ObjectToMap(STDFlowPath, MovedObjX, MovedObjY, FLOWgeoref, MovedMapX, MovedMapY);
Displacement3D(OrigMapX, OrigMapY, 0, MovedMapX, MovedMapY, 0, tempDist, azimuth, elevation);
if (tempDist > maxSnapLineDist) {
TableReadAttachment(SampleTBL, i, records, "point");
tempID1 = TableReadFieldStr(SampleTBL, sampleIDFldName, records[1]);
warnCountOutsideSearch = warnCountOutsideSearch + 1;
message = "Warning: Point " + tempID1 + " was found to be " + NumToStr(tempDist) + " meters from the closest computed standard\n" +
" flowpath line. This sample point will not be included in the remaining analysis\n" +
" of this script. Adjust the Flow Path and Basin parameters in the main dialog\n" +
" or modify the original location of the sample point and run the script again.\n\n";
print(message);
fwritestring(summaryFile, message);
} else {
# If user specifies, snap to the downstream node of the closet line within specified maximum search distance
snappedToNode = 0;
if (mainDLG.GetCtrlByID("id_SNAP_NODE_TOGGLE").GetValueNum()) {
P = GetVectorLine(STDFlowPath, lineNum);
# Get the two endpoints of the line and determine the distance to them from the orignal point
pt1 = P.GetVertex(0);
ObjectToMap(STDFlowPath, pt1.x, pt1.y, FLOWgeoref, Temp1MapX, Temp1MapY);
node1 = FindClosestNode(STDFlowPath, Temp1MapX, Temp1MapY, FLOWgeoref);
Displacement3D(MovedMapX, MovedMapY, 0, Temp1MapX, Temp1MapY, 0, pt1Dist, azimuth, elevation);
pt2 = P.GetVertex(P.GetNumPoints() - 1);
ObjectToMap(STDFlowPath, pt2.x, pt2.y, FLOWgeoref, Temp2MapX, Temp2MapY);
node2 = FindClosestNode(STDFlowPath, Temp2MapX, Temp2MapY, FLOWgeoref);
Displacement3D(MovedMapX, MovedMapY, 0, Temp2MapX, Temp2MapY, 0, pt2Dist, azimuth, elevation);
# If the downstream node is node1 and is within the snapping distance
if ( (pt1Dist < maxSnapNodeDist) and (node1 == STDFlowPath.line[lineNum].(LineToLineTblName).(downNodeFldName)) ) {
MovedMapX = Temp1MapX;
MovedMapY = Temp1MapY;
snappedToNode = node1;
# If the downstream node is node 2 and is within the snapping distance
} else if ( (pt2Dist < maxSnapNodeDist) and (node2 == STDFlowPath.line[lineNum].(LineToLineTblName).(downNodeFldName)) ) {
MovedMapX = Temp2MapX;
MovedMapY = Temp2MapY;
snappedToNode = node2;
}
}
# Snap the new points to the center of the DEM raster cells
MapToObject(DEMgeoref, MovedMapX, MovedMapY, SampleDEM, MovedObjX, MovedObjY);
MovedObjX = floor(MovedObjX) + 0.5;
MovedObjY = floor(MovedObjY) + 0.5;
ObjectToMap(SampleDEM, MovedObjX, MovedObjY, DEMgeoref, MovedMapX, MovedMapY);
# Convert from map coordinates to vector coordinates
MapToObject(WARPgeoref, MovedMapX, MovedMapY, WarpVECT, VectObjX, VectObjY);
# Check to see if the point already exists
pointExists = 0;
for j = 1 to numPointsAdded {
Temp1ObjX = ResultVECT.Point[j].Internal.x; # In object coordinates
Temp1ObjY = ResultVECT.Point[j].Internal.y;
# If both are in the center of the same cell, it's a duplicate
if ( (VectObjX == Temp1ObjX) and (VectObjY == Temp1ObjY) ) {
pointExists = j;
j = numPointsAdded;
}
}
# If the point already exists, mark the current point as a duplicate of the existing point
# The element number with the lowest value will always be the point added. Elements
# with a higher value will always be considered the duplicate since a point already exists.
if (pointExists > 0) {
TableReadAttachment(SampleTBL, i, records, "point");
tempID1 = TableReadFieldStr(SampleTBL, sampleIDFldName, records[1]);
TableReadAttachment(ResultSampleTBL, pointExists, records, "point");
tempID2 = TableReadFieldStr(ResultSampleTBL, sampleIDFldName, records[1]);
warnCountDuplicates = warnCountDuplicates + 1;
message = "Warning: Point " + tempID1 + " was found to be a duplicate of point " + tempID2 + " because both points\n" +
" were going to be placed in the center of the same cell of the DEM.\n" +
" The sample number has been recorded as being a duplicate of point " + tempID2 + ".\n\n";
print(message);
fwritestring(summaryFile, message);
TableReadAttachment(ResultPointToPointTBL, pointExists, records, "point");
currNumDups = TableReadFieldNum(ResultPointToPointTBL, numDuplFldName, records[1]);
currNumDups = currNumDups + 1;
# Add additional duplicate sample fields if additional are needed
while (currNumDups > currMaxNumDups) {
currMaxNumDups = currMaxNumDups + 1;
tempFieldName = dupSampleIDFldName + NumToStr(currMaxNumDups);
# Insert the field before the "OnLine" field
TableInsertFieldString(ResultPointToPointTBL, tempFieldName, beforeField, 12);
}
TableReadAttachment(SampleTBL, i, records, "point");
sampleNum$ = TableReadFieldStr(SampleTBL, sampleIDFldName, records[1]);
# Add the sample point to the table as a duplicate
TableReadAttachment(ResultPointToPointTBL, pointExists, records, "point");
TableWriteField(ResultPointToPointTBL, records[1], numDuplFldName, currNumDups);
tempFieldName = dupSampleIDFldName + NumToStr(currNumDups);
TableWriteField(ResultPointToPointTBL, records[1], tempFieldName, sampleNum$);
# Add the point if it doesn't exist
} else {
numPointsAdded = numPointsAdded + 1;
x[numPointsAdded] = MovedObjX - 1;
y[numPointsAdded] = MovedObjY - 1;
# Check if the point ended up getting moved to a downstream node. If so, attach it to that downstream line
node1 = FindClosestNode(STDFlowPath, MovedMapX, MovedMapY, FLOWgeoref);
pt1.x = STDFlowPath.node[node1].Internal.x;
pt1.y = STDFlowPath.node[node1].Internal.y;
pt1 = ObjectToMap(STDFlowPath, pt1.x, pt1.y, FLOWgeoref);
pt1 = MapToObject(DEMgeoref, pt1.x, pt1.y, SampleDEM);
pt1.x = floor(pt1.x) + 0.5;
pt1.y = floor(pt1.y) + 0.5;
# If the node and the point are in the center of the same raster cell, then consider the point on the downstream
# line of that node
if ((pt1.x == MovedObjX) && (pt1.y == MovedObjY)) {
# If the closest node is the downstream node
if (node1 == STDFlowPath.line[lineNum].(LineToLineTblName).(downNodeFldName)) {
# Assign the line the sample is on to the downstream line
if (STDFlowPath.line[lineNum].(LineToLineTblName).(downLineFldName) != -1) {
lineNum = STDFlowPath.line[lineNum].(LineToLineTblName).(downLineFldName);
}
}
}
# Assumes the vector point added has the element number of the current numPointsAdded
VectorAddPoint(ResultVECT, VectObjX, VectObjY, SampleDEM[MovedObjX, MovedObjY]);
numRecords = TableReadAttachment(SampleTBL, i, records, "point");
err = TableWriteAttachment(ResultSampleTBL, numPointsAdded, records, numRecords, "point");
if (err <= 0) {
message = "Warning: Either an error occured while attaching records to the new point vector,\n" +
" or there are no attached records for element point number: " + NumToStr(i) + ".\n\n";
print(message);
fwritestring(summaryFile, message);
}
# Calculate the distance the sample point was displaced by
Displacement3D(OrigMapX, OrigMapY, 0, MovedMapX, MovedMapY, 0, distance, azimuth, elevation);
# Add the moved points and displacement values to the database
for j = 1 to numRecords {
TableWriteField(ResultSampleTBL, records[j], modEastFldName, MovedMapX);
TableWriteField(ResultSampleTBL, records[j], modNorthFldName, MovedMapY);
TableWriteField(ResultSampleTBL, records[j], displacedFldName, distance);
TableWriteField(ResultSampleTBL, records[j], snapToDownNodeFldName, snappedToNode);
}
records[1] = TableWriteRecord(ResultPointToPointTBL, 0, lineNum);
TableWriteAttachment(ResultPointToPointTBL, numPointsAdded, records, 1, "point");
}
}
}
VectorValidate(ResultVECT);
if (NumVectorPoints(ResultVECT) == 0) {
message = "None of the sample points were found to be within the specified maximum search distance.\n" +
"Please verify that your search distance settings are appropriate.\n \n" +
"Hint: The standard flowpath vector has already been created and should be in the output project file.\n" +
" Use this flowpath vector object along with your point vector to determine appropriate setting values.";
PopupMessage(message);
print(message);
fwritestring(summaryFile, message);
CloseVector(ResultVECT);
CloseRaster(SampleDEM);
CloseVector(SampleVECT);
CloseVector(WarpVECT);
return;
}
message = "##########################################################################################\n" +
"Number of Unmovable points found: " + NumToStr(warnCountOutsideSearch) + "\n" +
"Number of Duplicate points found: " + NumToStr(warnCountDuplicates) + "\n" +
"##########################################################################################\n\n";
print(message);
fwritestring(summaryFile, message);
#######################################################################################################
# Copy the multiple stream order schemas to the point sediment sample table
statusC.Message = "Step 6 of 18: Copying stream order schemas to sediment databases";
#######################################################################################################
numSamplePoints = NumVectorPoints(ResultVECT);
statusC.BarInit(numSamplePoints, 0);
for i = 1 to numSamplePoints {
statusC.BarIncrement(1, 0);
TableReadAttachment(ResultPointToPointTBL, i, records, "point");
lineNum = TableReadFieldNum(ResultPointToPointTBL, onLineFldName, records[1]);
TableReadAttachment(ResultSampleTBL, i, records, "point");
TableWriteField(ResultSampleTBL, records[1], hortonFldName, STDFlowPath.line[lineNum].STREAM_ORDER.Horton);
TableWriteField(ResultSampleTBL, records[1], shreveFldName, STDFlowPath.line[lineNum].STREAM_ORDER.Shreve);
TableWriteField(ResultSampleTBL, records[1], strahlerFldName, STDFlowPath.line[lineNum].STREAM_ORDER.Strahler);
TableWriteField(ResultSampleTBL, records[1], scheideggerFldName, STDFlowPath.line[lineNum].STREAM_ORDER.Scheidegger);
}
#######################################################################################################
# Caluculate the upstream channel length for each sample point
statusC.Message = "Step 7 of 18: Computing upstream channel length for each sample point";
#######################################################################################################
local numeric closestVertex, testVertex;
local numeric distCurr2Test, distClose2Test, distCurr2Close, tempDist;
local numeric lengthClose2Node, lengthCurr2Node, upstreamChanLength, tempLength, junk1, junk2;
local class POINT2D currPt, tempPt1, tempPt2, testPt, closePt;
statusC.BarInit(numSamplePoints, 0);
for i = 1 to numSamplePoints {
statusC.BarIncrement(1, 0);
lengthCurr2Node = 0;
lengthClose2Node = 0;
tempLength = 0;
TableReadAttachment(ResultPointToPointTBL, i, records, "point");
lineNum = TableReadFieldNum(ResultPointToPointTBL, onLineFldName, records[1]);
P = GetVectorLine(STDFlowPath, lineNum);
# Get one of the endpoints
pt1 = P.GetVertex(0);
ObjectToMap(STDFlowPath, pt1.x, pt1.y, FLOWgeoref, Temp1MapX, Temp1MapY);
node1 = FindClosestNode(STDFlowPath, Temp1MapX, Temp1MapY, FLOWgeoref);
# Node 1 needs to be upstream. If it's not, turn the Polyline around
if (node1 != STDFlowPath.line[lineNum].(LineToLineTblName).(upNodeFldName)) {
P.Reverse();
pt1 = P.GetVertex(0);
}
pt1 = ObjectToMap(STDFlowPath, pt1.x, pt1.y, FLOWgeoref);
# Find the closest vertex on the line to the current sample point
TableReadAttachment(ResultSampleTBL, i, records, "point");
currPt.x = TableReadFieldNum(ResultSampleTBL, modEastFldName, records[1]);
currPt.y = TableReadFieldNum(ResultSampleTBL, modNorthFldName, records[1]);
currPt = MapToObject(FLOWgeoref, currPt.x, currPt.y, STDFlowPath);
closestVertex = P.FindClosestVertex(currPt);
currPt = ObjectToMap(STDFlowPath, currPt.x, currPt.y, FLOWgeoref);
# If closestVertex is the upstream node, then the line distance is the distance from the current
# point to the upstream node;
if (closestVertex == 0) {
Displacement3D(currPt.x, currPt.y, 0, pt1.x, pt1.y, 0, lengthCurr2Node, junk1, junk2);
# Otherwise, need to find distance from closest vertex to the upnode +/- where the actual point is
} else {
# Calculate the length of the line from the closest vertex to the upstream node
for j = 1 to closestVertex {
tempPt1 = P.GetVertex(j-1);
tempPt1 = ObjectToMap(STDFlowPath, tempPt1.x, tempPt1.y, FLOWgeoref);
tempPt2 = P.GetVertex(j);
tempPt2 = ObjectToMap(STDFlowPath, tempPt2.x, tempPt2.y, FLOWgeoref);
Displacement3D(tempPt1.x, tempPt1.y, 0, tempPt2.x, tempPt2.y, 0, tempDist, junk1, junk2);
lengthClose2Node = lengthClose2Node + tempDist;
}
closePt = P.GetVertex(closestVertex);
closePt = ObjectToMap(STDFlowPath, closePt.x, closePt.y, FLOWgeoref);
# If the closest vertex point is not the same as the current point, offset the vertex to node
# length based on the distance the current point is from the closest vertex
if (closePt != currPt) {
# Get the vertex that is one closer to the upstream node
testVertex = closestVertex - 1;
testPt = P.GetVertex(testVertex);
testPt = ObjectToMap(STDFlowPath, testPt.x, testPt.y, FLOWgeoref);
Displacement3D(currPt.x, currPt.y, 0, testPt.x, testPt.y, 0, distCurr2Test, junk1, junk2);
Displacement3D(closePt.x, closePt.y, 0, testPt.x, testPt.y, 0, distClose2Test, junk1, junk2);
Displacement3D(currPt.x, currPt.y, 0, closePt.x, closePt.y, 0, distCurr2Close, junk1, junk2);
# If the distance from the current point to the test vertex is greater then the closest to the test,
# then the additional distance needs to be added.
if (distCurr2Test > distClose2Test) {
lengthCurr2Node = lengthClose2Node + distCurr2Close;
# Otherwise, the distance needs to be subtracted
} else {
lengthCurr2Node = lengthClose2Node - distCurr2Close;
}
} else {
lengthCurr2Node = lengthClose2Node;
}
}
# Get the total upstream channel length that is upstream of this line
TableReadAttachment(STDFlowPathTBL, lineNum, records, "line");
upstreamChanLength = TableReadFieldNum(STDFlowPathTBL, totalUpChanLengthFldName, records[1]);
# Compute the total upstream channel length from the sample point in km
tempLength = upstreamChanLength + (lengthCurr2Node / 1000);
# Write the length to the record
TableReadAttachment(ResultSampleTBL, i, records, "point");
TableWriteField(ResultSampleTBL, records[1], totalUpChanLenFldName, tempLength);
} # End for i loop
#######################################################################################################
# Create relationships between sample points that are on the same line based upon their location
# upstream and downstream of each other
statusC.Message = "Step 8 of 18: Computing relationships between sample points on same flowpath line";
#######################################################################################################
local numeric currVertex, tempVertex;
local numeric closestSample, closestDist, distApart;
local numeric numOnSameLine, thisLineNum;
local numeric tempDist, currDist, closeDist;
local array numeric onSameLine[10]; # Currently allows up to 10 points on the same line
# Create relationships for points that ARE on the same line as others
statusC.BarInit(numSamplePoints, 0);
for i = 1 to numSamplePoints {
statusC.BarIncrement(1, 0);
thisLineNum = ResultVECT.point[i].PointToPoint.(onLineFldName);
closestDist = -1;
closestSample = -1;
closestVertex = -1;
# Search for any points that are on the same line
for j = 1 to numSamplePoints {
# Reset variables
numOnSameLine = 0;
for k = 1 to 10 {
onSameLine[k] = -1;
}
# Ignore itself
if (j != i) {
# Check if the other point is on the same line
if (thisLineNum == ResultVECT.point[j].PointToPoint.(onLineFldName)) {
numOnSameLine = numOnSameLine + 1;
onSameLine[numOnSameLine] = j;
}
}
# If there are points on the same line
if (numOnSameLine > 0) {
# Get one of the end nodes and check if it is upstream
P = GetVectorLine(STDFlowPath, thisLineNum);
# Get one of the endpoints
pt1 = P.GetVertex(0);
ObjectToMap(STDFlowPath, pt1.x, pt1.y, FLOWgeoref, Temp1MapX, Temp1MapY);
node1 = FindClosestNode(STDFlowPath, Temp1MapX, Temp1MapY, FLOWgeoref);
# Node 1 needs to be downstream. If it's not, turn the Polyline around
if (node1 != STDFlowPath.line[thisLineNum].(LineToLineTblName).(downNodeFldName)) {
P.Reverse();
pt1 = P.GetVertex(0);
}
# Find the closest vertex on the line to the current sample point
TableReadAttachment(ResultSampleTBL, i, records, "point");
currPt.x = TableReadFieldNum(ResultSampleTBL, modEastFldName, records[1]);
currPt.y = TableReadFieldNum(ResultSampleTBL, modNorthFldName, records[1]);
currPt = MapToObject(FLOWgeoref, currPt.x, currPt.y, STDFlowPath);
currVertex = P.FindClosestVertex(currPt);
for k = 1 to numOnSameLine {
TableReadAttachment(ResultSampleTBL, onSameLine[k], records, "point");
tempPt1.x = TableReadFieldNum(ResultSampleTBL, modEastFldName, records[1]);
tempPt1.y = TableReadFieldNum(ResultSampleTBL, modNorthFldName, records[1]);
tempPt1 = MapToObject(FLOWgeoref, tempPt1.x, tempPt1.y, STDFlowPath);
tempVertex = P.FindClosestVertex(tempPt1);
# Notify user if sample points are located at the same place. Shouldn't happen.
if (tempPt1 == currPt) {
TableReadAttachment(ResultSampleTBL, onSameLine[k], records, "point");
tempID1 = TableReadFieldStr(ResultSampleTBL, sampleIDFldName, records[1]);
TableReadAttachment(ResultSampleTBL, i, records, "point");
tempID2 = TableReadFieldStr(ResultSampleTBL, sampleIDFldName, records[1]);
message = "Warning: Sample point " + tempID1 + " has the same location as " + tempID2 + " \n" +
" It is best that these points be edited and then run again.\n\n";
print(message);
fwritestring(summaryFile, message);
# Check if tempVertex could be closer then closestVertex and is the same or less then the current
} else if ( (tempVertex >= closestVertex) and (tempVertex <= currVertex) ) {
# If the temp Vertex is the same as the current Vertex, check if temp is downstream
if (tempVertex == currVertex) {
# Use the next vertex downstream for the test
testVertex = currVertex - 1;
# Make sure the next vertex downstream is valid. Can't be less then 0 or the
# point would be on the next downstream line.
if (testVertex < 0) {
testVertex = 0;
}
# Determine the distance from each point and the distances to the next downtream node
testPt = P.GetVertex(testVertex);
Displacement3D(currPt.x, currPt.y, 0, testPt.x, testPt.y, 0, currDist, junk1, junk2);
Displacement3D(tempPt1.x, tempPt1.y, 0, testPt.x, testPt.y, 0, tempDist, junk1, junk2);
Displacement3D(currPt.x, currPt.y, 0, tempPt1.x, tempPt1.y, 0, distApart, junk1, junk2);
# If the temp point is closer to the downstream vertex
if (tempDist < currDist) {
# If this is the first temp at this vertex
if (closestDist == -1) {
closestDist = distApart;
closestSample = onSameLine[k];
closestVertex = tempVertex;
# If another point downstream had the same vertex, but this one is closer
} else if ( (closestDist != -1) and (distApart < closestDist) ) {
closestDist = distApart;
closestSample = onSameLine[k];
closestVertex = tempVertex;
}
}
# If the temp Vertex is downstream of the current Vertex
} else if (tempVertex < currVertex) {
# If the temp Vertex is the same as the closest vertex, find out which is further upstream
if (tempVertex == closestVertex) {
# Use the next vertex upstream as a test
testVertex = closestVertex + 1;
# Get the current closests points coordinates
TableReadAttachment(ResultSampleTBL, closestSample, records, "point");
closePt.x = TableReadFieldNum(ResultSampleTBL, modEastFldName, records[1]);
closePt.y = TableReadFieldNum(ResultSampleTBL, modNorthFldName, records[1]);
closePt = MapToObject(FLOWgeoref, closePt.x, closePt.y, STDFlowPath);
# Determine the distance from each point and the distances to the next uptream node
testPt = P.GetVertex(testVertex);
Displacement3D(closePt.x, closePt.y, 0, testPt.x, testPt.y, 0, closeDist, junk1, junk2);
Displacement3D(tempPt1.x, tempPt1.y, 0, testPt.x, testPt.y, 0, tempDist, junk1, junk2);
Displacement3D(closePt.x, closePt.y, 0, tempPt1.x, tempPt1.y, 0, distApart, junk1, junk2);
# If the temp point is closer to the upstream vertex, then it is further upstream
if (tempDist < closeDist) {
closestSample = onSameLine[k];
closestVertex = tempVertex;
}
# If temp Vertex is greater then closest Vertex, then it is closer to the current point
} else if (tempVertex > closestVertex) {
closestSample = onSameLine[k];
closestVertex = tempVertex;
}
} # End (temp vs curr)
} # End (temp >= closest)
} # End for loop
# Recored the sample point number found upstream
if (closestVertex != -1) {
TableReadAttachment(ResultSampleTBL, closestSample, records, "point");
sampleNum$ = TableReadFieldStr(ResultSampleTBL, sampleIDFldName, records[1]);
TableReadAttachment(ResultPointToPointTBL, i, records, "point");
TableWriteField(ResultPointToPointTBL, records[1], downSampleFldName, sampleNum$);
}
} # End "For each on same line"
} #End "On Same Line"
}
#######################################################################################################
# Determine the most upstream and downstream sample points on each line segment
statusC.Message = "Step 9 of 18: Computing the most upstream and downstream sample points for each flowpath line";
#######################################################################################################
local array numeric pointsOnLine[100];
local numeric numPointsOnLine, relationshipExists;
local string currSampleNum$;
# Determine the most upstream and most downstream sample points on each line
statusC.BarInit(numFlowPathLines, 0);
for i = 1 to numFlowPathLines {
statusC.BarIncrement(1, 0);
numPointsOnLine = 0;
# For each sample point
for j = 1 to numSamplePoints {
# If the sample point is on the current line, add it to the array
if (ResultVECT.point[j].PointToPoint.(onLineFldName) == i) {
numPointsOnLine = numPointsOnLine + 1;
pointsOnLine[numPointsOnLine] = j;
}
}
# If the line has no sample points on it, leave the fields blank
if (numPointsOnLine == 0) {
# Assign it a value of nothing, can be changed to a unique identifier string
TableReadAttachment(STDFlowPathTBL, i, records, "line");
TableWriteField(STDFlowPathTBL, records[1], mostDownSampleFldName, "");
TableWriteField(STDFlowPathTBL, records[1], mostUpSampleFldName, "");
# If there is only one point on the line, it is both the most upstream and downstream sample
} else if (numPointsOnLine == 1) {
TableReadAttachment(ResultSampleTBL, pointsOnLine[1], records, "point");
sampleNum$ = TableReadFieldStr(ResultSampleTBL, sampleIDFldName, records[1]);
TableReadAttachment(STDFlowPathTBL, i, records, "line");
TableWriteField(STDFlowPathTBL, records[1], mostDownSampleFldName, sampleNum$);
TableWriteField(STDFlowPathTBL, records[1], mostUpSampleFldName, sampleNum$);
# If there are more then 1 sample points on the line
} else if (numPointsOnLine > 1) {
# For each point on the line, determine it's most downstream sample
for j = 1 to numPointsOnLine {
# If the current point on the line doesn't have a Downstream Sample, it is the most downstream
TableReadAttachment(ResultPointToPointTBL, pointsOnLine[j], records, "point");
sampleNum$ = TableReadFieldStr(ResultPointToPointTBL, downSampleFldName, records[1]);
if (sampleNum$.length == 0) {
# Get it's real sample number and set it as the most downstream sample for this line
TableReadAttachment(ResultSampleTBL, pointsOnLine[j], records, "point");
sampleNum$ = TableReadFieldStr(ResultSampleTBL, sampleIDFldName, records[1]);
TableReadAttachment(STDFlowPathTBL, i, records, "line");
TableWriteField(STDFlowPathTBL, records[1], mostDownSampleFldName, sampleNum$);
}
}
# For each point on the line, determine it's most upstream sample
for j = 1 to numPointsOnLine {
# If no points relate to it as being downstream, it is the most upstream in the line
relationshipExists = 0;
TableReadAttachment(ResultSampleTBL, pointsOnLine[j], records, "point");
currSampleNum$ = TableReadFieldStr(ResultSampleTBL, sampleIDFldName, records[1]);
# Check the other points on the line to see if the point to it as being downstream
for k = 1 to numPointsOnLine {
# Don't check itself
if (k != j) {
TableReadAttachment(ResultPointToPointTBL, pointsOnLine[k], records, "point");
sampleNum$ = TableReadFieldStr(ResultPointToPointTBL, downSampleFldName, records[1]);
# If something is pointing to it as downstream, it isn't the most upstream
if (sampleNum$ == currSampleNum$) {
relationshipExists = 1;
}
}
}
# If no relationship exists on this line, it is the most upstream point on that line
if (!relationshipExists) {
TableReadAttachment(ResultSampleTBL, pointsOnLine[j], records, "point");
sampleNum$ = TableReadFieldStr(ResultSampleTBL, sampleIDFldName, records[1]);
TableReadAttachment(STDFlowPathTBL, i, records, "line");
TableWriteField(STDFlowPathTBL, records[1], mostUpSampleFldName, sampleNum$);
}
} #End For loop
} # End comparision of numPointsOnLine
} # End for loop
#######################################################################################################
# Determine the relationships for sample points that require traversing a seperate line to find their
# immediate downstream or multiple upstream sample points
statusC.Message = "Step 10 of 18: Searching for downstream sample points for each sample";
#######################################################################################################
local numeric currLineNum, downstreamLineNum, stillSearching;
# Find the downstream relationships first
statusC.BarInit(numSamplePoints, 0);
for i = 1 to numSamplePoints {
statusC.BarIncrement(1, 0);
TableReadAttachment(ResultPointToPointTBL, i, records, "point");
sampleNum$ = TableReadFieldStr(ResultPointToPointTBL, downSampleFldName, records[1]);
# If the sample doesn't currently have a downstream sample, it needs to find one
if (sampleNum$.length == 0) {
TableReadAttachment(ResultPointToPointTBL, i, records, "point");
currLineNum = TableReadFieldNum(ResultPointToPointTBL, onLineFldName, records[1]);
TableReadAttachment(STDFlowPathTBL, currLineNum, records, "line");
downstreamLineNum = TableReadFieldNum(STDFlowPathTBL, downLineFldName, records[1]);
# Check if the current line has a downstream line. If so, traverse downstream until we find a
# line with an UpstreamSample, if it exists
if (downstreamLineNum != -1) {
stillSearching = 1;
while (stillSearching) {
TableReadAttachment(STDFlowPathTBL, downstreamLineNum, records, "line");
sampleNum$ = TableReadFieldStr(STDFlowPathTBL, mostUpSampleFldName, records[1]);
# If the downstream line has a Upstream Sample, attach it as the downstream sample
if (sampleNum$.length != 0) {
# Write the record and stop searching
TableReadAttachment(ResultPointToPointTBL, i, records, "point");
TableWriteField(ResultPointToPointTBL, records[1], downSampleFldName, sampleNum$);
stillSearching = 0;
# Otherwise, check to see if the next downstream line exists
} else {
currLineNum = downstreamLineNum;
TableReadAttachment(STDFlowPathTBL, currLineNum, records, "line");
downstreamLineNum = TableReadFieldNum(STDFlowPathTBL, downLineFldName, records[1]);
# If it doesn't exist, stop searching for it
if (downstreamLineNum == -1) {
stillSearching = 0;
}
# If it does exist, repeat the while loop to search that downstream line
} # End sampleNum$.length
} # End while loop
} # End initial downstream test
} # End downstream sample check
} # End for each point loop
local numeric currMaxNumUpstreamSampleFields = 2;
local array numeric upstreamSamples[numSamplePoints];
local numeric numUpstreamSamples;
# Find the upstream relationships between sample points
statusC.BarInit(numSamplePoints, 0);
for i = 1 to numSamplePoints {
statusC.BarIncrement(1, 0);
numUpstreamSamples = 0;
TableReadAttachment(ResultSampleTBL, i, records, "point");
currSampleNum$ = TableReadFieldStr(ResultSampleTBL, sampleIDFldName, records[1]);
# Check all samples to see if they list it as downstream of them
for j = 1 to numSamplePoints {
# Don't bother check itself
if (j != i) {
TableReadAttachment(ResultPointToPointTBL, j, records, "point");
sampleNum$ = TableReadFieldStr(ResultPointToPointTBL, downSampleFldName, records[1]);
# If the j sample point's downstream sample is i, then save j to the array
if (sampleNum$ == currSampleNum$) {
numUpstreamSamples = numUpstreamSamples + 1;
upstreamSamples[numUpstreamSamples] = j;
} # End sample comparison
} # End if not itself
} # End for j loop
# If there aren't enough upstream sample fields in the table, create the additional number needed
while (currMaxNumUpstreamSampleFields < numUpstreamSamples) {
currMaxNumUpstreamSampleFields = currMaxNumUpstreamSampleFields + 1;
tempFieldName = upSampleFldName + NumToStr(currMaxNumUpstreamSampleFields);
TableAddFieldString(ResultPointToPointTBL, tempFieldName, 11);
}
TableReadAttachment(ResultPointToPointTBL, i, records, "point");
TableWriteField(ResultPointToPointTBL, records[1], numUpSamplesFldName, numUpstreamSamples);
# Add each upstream sample number to the current samples database table
for j = 1 to numUpstreamSamples {
# Get the upstream sample number
TableReadAttachment(ResultSampleTBL, upstreamSamples[j], records, "point");
sampleNum$ = TableReadFieldStr(ResultSampleTBL, sampleIDFldName, records[1]);
# Add the sample number as being upstream in the appropriate field
TableReadAttachment(ResultPointToPointTBL, i, records, "point");
tempFieldName = upSampleFldName + NumToStr(j);
TableWriteField(ResultPointToPointTBL, records[1], tempFieldName, sampleNum$);
} # End for j loop
} # End for i loop
#######################################################################################################
# Find the total number of upstream sample points for each sample point
statusC.Message = "Step 11 of 18: Computing total number of upstream sample points for each sample";
#######################################################################################################
local numeric totalUpstream, numUpstream, upstreamSampleTotal, currTotalUpSamples;
numComputed = 0;
tempCount = 0;
# Initialize the total number of upstream samples field by setting it to 0 if there are non upstream,
# or to -1 to specify that it needs to be computed.
statusC.BarInit(numSamplePoints, 0);
for i = 1 to numSamplePoints {
statusC.BarIncrement(1, 0);
TableReadAttachment(ResultPointToPointTBL, i, records, "point");
totalUpstream = TableReadFieldNum(ResultPointToPointTBL, numUpSamplesFldName, records[1]);
if (totalUpstream == 0) {
numComputed = numComputed + 1;
TableWriteField(ResultPointToPointTBL, records[1], totalUpSamplesFldName, 0);
} else {
TableWriteField(ResultPointToPointTBL, records[1], totalUpSamplesFldName, -1);
}
}
# Cycle through all the points continually until they all have the correct number of total upstream samples
statusC.BarInit(numSamplePoints, 0);
while (numComputed > 0 ) {
# Reset the while counter
tempCount += numComputed;
statusC.BarIncrement(numComputed, 0);
numComputed = 0;
for i = 1 to numSamplePoints {
# Check the number of total upstream samples
TableReadAttachment(ResultPointToPointTBL, i, records, "point");
totalUpstream = TableReadFieldNum(ResultPointToPointTBL, totalUpSamplesFldName, records[1]);
# If -1, it still needs to be computed
if (totalUpstream == -1) {
# Find the number of immediate upstream samples and add it to the current total
numUpstream = TableReadFieldNum(ResultPointToPointTBL, numUpSamplesFldName, records[1]);
currTotalUpSamples = numUpstream;
# Initialize the list and then get the list of the upstream sample numbers
UpSample.Clear();
UpSample.AddToFront("");
for j = 1 to numUpstream {
UpSample.AddToEnd(TableReadFieldStr(ResultPointToPointTBL, upSampleFldName + NumToStr(j), records[1]));
}
for j = 1 to numUpstream {
# If the curret Total is not invalid, keep checking upstream samples
if (currTotalUpSamples != -1) {
# Get the point element the Sample Number refers to and check it's total upstream samples field
records[1] = TableKeyFieldLookup(ResultSampleTBL, sampleIDFldName, UpSample[j], "Equals");
TableGetRecordElementList(ResultSampleTBL, records[1], elements, "point");
TableReadAttachment(ResultPointToPointTBL, elements[1], records, "point");
upstreamSampleTotal = TableReadFieldNum(ResultPointToPointTBL, totalUpSamplesFldName, records[1]);
# If the upstream sample point of i has a total number of sample points, add it to the running total
if (upstreamSampleTotal != -1) {
currTotalUpSamples = currTotalUpSamples + upstreamSampleTotal;
# The upstream samples have yet to be calculated, set the current total to invalid
} else {
currTotalUpSamples = -1;
}
}
} # End for j loop
# If a valid current total of upstream samples was collected, record it to the table for that point
if (currTotalUpSamples != -1) {
TableReadAttachment(ResultPointToPointTBL, i, records, "point");
TableWriteField(ResultPointToPointTBL, records[1], totalUpSamplesFldName, currTotalUpSamples);
# Update the while counter so it will run again
numComputed = numComputed + 1;
}
} # End if
} # End for i loop
} # End while loop
# Record the number of sample upstream to the point sediment sample database
statusC.BarInit(numSamplePoints, 0);
for i = 1 to numSamplePoints {
statusC.BarIncrement(1, 0);
TableReadAttachment(ResultPointToPointTBL, i, records, "point");
totalUpstream = TableReadFieldNum(ResultPointToPointTBL, totalUpSamplesFldName, records[1]);
TableReadAttachment(ResultSampleTBL, i, records, "point");
TableWriteField(ResultSampleTBL, records[1], totalUpSamplesFldName, totalUpstream);
}
#######################################################################################################
# Use the new sample point locations to compute a new user defined watershed
statusC.Message = "Step 12 of 18 Computing user defined Watershed objects";
statusC.BarClear(0);
#######################################################################################################
# Recreate the watershed using the new modified points
string watershedparms$ = "";
if ( mainDLG.GetCtrlByID("id_FILL_ALL_DEP").GetValueNum()) {
watershedparms$ += "FillAllDepressions,";
}
if ( mainDLG.GetCtrlByID("id_COMP_USR_FLOWPATH").GetValueNum()) {
watershedparms$ += "FlowPath,";
}
# if ( mainDLG.GetCtrlByID("id_COMP_STD_DBASE").GetValueNum()) {
# watershedparms$ += "Database,";
# }
watershedparms$ += "Basin";
WatershedComputeElements(watershed, x, y, numSamplePoints, watershedparms$);
if ( mainDLG.GetCtrlByID("id_COMP_USR_FLOWPATH").GetValueNum()) {
WatershedGetObject(watershed, "VectorUserFlowPath", ObjectFilename, ObjectName);
ObjNumber = ObjectNumber(ObjectFilename, ObjectName, "VECTOR");
CopyObject(ObjectFilename, ObjNumber, OutputProjectFile);
}
WatershedGetObject(watershed, "VectorUserBasin", ObjectFilename, ObjectName);
ObjNumber = ObjectNumber(ObjectFilename, ObjectName, "VECTOR");
CopyObject(ObjectFilename, ObjNumber, OutputProjectFile);
# Open the "VectorUserBasin" and get it's georef
OpenVector(UBasinVECT, OutputProjectFile, "USDBASIN");
local class GEOREF UBSNgeoref = GetLastUsedGeorefObject(UBasinVECT);
#######################################################################################################
# Attach the user basins with the appropriate sample points that are within the basin
statusC.Message = "Step 13 of 18: Attaching catchment basin polygons to appropriate sample points";
statusC.BarClear(0);
#######################################################################################################
local numeric numBasinsComputed = 0;
local numeric numUserBasins = NumVectorPolys(UBasinVECT);
local numeric polyNum;
local class DATABASE UBasinDB = OpenVectorPolyDatabase(UBasinVECT);
# Copy the sample database table from original vector
if (TableCopy(SampleDB, SampleTBL, UBasinDB) < 0) {
message = "Error occured while trying to copy the sample table.\n" + "Exiting the script now... \n";
print(message);
fwritestring(summaryFile, message);
CloseVector(UBasinVECT);
CloseRaster(SampleDEM);
CloseVector(SampleVECT);
CloseVector(WarpVECT);
CloseVector(ResultVECT);
return;
}
local class DBTABLEINFO UBasinSampleTBL = DatabaseGetTableInfo(UBasinDB, sampleTblName);
UBasinSampleTBL.OneRecordPerElement = 1;
local class DBTABLEINFO UBasinPolyToPolyTBL = TableCreate(UBasinDB, PolyToPolyTblName, PolyToPolyTblDesc);
UBasinPolyToPolyTBL.OneRecordPerElement = 1;
TableAddFieldString(UBasinPolyToPolyTBL, sampleIDFldName, 12); # Sample ID realated to the basin polygon
TableAddFieldInteger(UBasinPolyToPolyTBL, numPolysPerSampleFldName, 15); # Number of vector polygons related to that sample
TableAddFieldFloat(UBasinPolyToPolyTBL, areaSamplePolysFldName, 15, 4); # Total of all vector polygons related to that sample
# Create a record in the table for each basin
statusC.BarInit(numUserBasins, 0);
for i = 1 to numUserBasins {
statusC.BarIncrement(1, 0);
records[1] = TableWriteRecord(UBasinPolyToPolyTBL, 0, "", 1, -1);
TableWriteAttachment(UBasinPolyToPolyTBL, i, records, 1, "polygon");
if (DEBUG) {
fwritestring(summaryFile, " i:" + NumToStr(i) + " rec[1]: " + NumToStr(records[1]) + "\n");
}
}
statusC.BarInit(numSamplePoints, 0);
for i = 1 to numSamplePoints {
statusC.BarIncrement(1, 0);
Temp1ObjX = ResultVECT.Point[i].Internal.x; # In object coordinates
Temp1ObjY = ResultVECT.Point[i].Internal.y;
ObjectToMap(ResultVECT, Temp1ObjX, Temp1ObjY, RSLTgeoref, Temp1MapX, Temp1MapY);
polyNum = FindClosestPoly(UBasinVECT, Temp1MapX, Temp1MapY, UBSNgeoref);
# If a polygon was found, save the sample name to the polygon table
if (polyNum > 0) {
TableReadAttachment(ResultSampleTBL, i, records, "point");
sampleNum$ = TableReadFieldStr(ResultSampleTBL, sampleIDFldName, records[1]);
TableReadAttachment(UBasinPolyToPolyTBL, polyNum, records, "polygon");
TableWriteField(UBasinPolyToPolyTBL, records[1], sampleIDFldName, sampleNum$);
TableWriteField(UBasinPolyToPolyTBL, records[1], numPolysPerSampleFldName, 1);
# Warning if no polygon is found. This shouldn't happen.
} else {
TableReadAttachment(ResultSampleTBL, i, records, "point");
tempID1 = TableReadFieldStr(ResultSampleTBL, sampleIDFldName, records[1]);
message = "Warning: The sample point with the element number " + NumToStr(i) + " was not able to be associated\n" +
" with one of the polygons. Please examine the results to confirm that the\n" +
" sample point does not fall within any of the user defined basin polygons.\n\n";
print(message);
fwritestring(summaryFile, message);
}
}
#######################################################################################################
# Remove records for basins that are not attached to any sample point, and warn the user of them.
statusC.Message = "Step 14 of 18: Removing unattached polygon records";
#######################################################################################################
numRecords = NumRecords(UBasinPolyToPolyTBL);
# Scan through the records. If it doesn't have a sample number, it shouldn't be a watershed basin.
# (Some extra polygons are created due to multiple polygons created around it)
statusC.BarInit(numRecords, 0);
for i = 1 to NumRecords(UBasinPolyToPolyTBL) {
statusC.BarIncrement(1, 0);
# Check if it has a sample number attached
sampleNum$ = TableReadFieldStr(UBasinPolyToPolyTBL, sampleIDFldName, i);
if (DEBUG) {
fwritestring(summaryFile, " i:" + NumToStr(i) + " samp$:" + sampleNum$ + "\n");
}
# If it doesn't, get the polygon element and remove it's attachment to the record
if (sampleNum$.length == 0) {
numeric tempnum;
tempnum = TableGetRecordElementList(UBasinPolyToPolyTBL, i, elements, "polygon");
records[1] = i;
TableRemoveAttachment(UBasinPolyToPolyTBL, elements[1], records, 1, "polygon");
if (DEBUG) {
fwritestring(summaryFile, " rec[1]:" + NumToStr(records[1]) + " tempN:" + NumToStr(tempnum) + " elem[1]:" + NumToStr(elements[1]) + "\n");
}
warnCountLostPolygons = warnCountLostPolygons + 1;
message = "Warning: The following polygon element " + NumToStr(elements[1]) + " does not have a sample point attached\n" +
" to it. This could be a result of a squeeze point in the watershed basins\n" +
" vector. Modification of the original point may be required.\n\n";
print(message);
fwritestring(summaryFile, message);
}
}
# Remove all the unattached records
TableRemoveUnattachedRecords(UBasinPolyToPolyTBL);
message = "##########################################################################################\n" +
"Number of unidentified polygons found: " + NumToStr(warnCountLostPolygons) + "\n" +
"##########################################################################################\n\n";
print(message);
fwritestring(summaryFile, message);
#######################################################################################################
# Search through the unattached basins and if they are connected to another basin, which has a SampleID,
# by only a single node and a flowpath line that intersects that node and the two polygons, assign it
# the same sample ID. This solves the case in which a catchment basin is "pinched" due to a diagonal
# flowpath line and DEM values that cause two or more seperate polygons to be created for that seed point.
# A special case which this will not find is if the "should be attached" basin was not large enough for
# a flowpath line to be created. Without the flowpath line, the relationship could not be made.
statusC.Message = "Step 15 of 18: Searching for relevant unattached polygons";
#######################################################################################################
local numeric m, n, t, u, tempValx, tempValy, stillPassingTest, numPolys;
local numeric numNodes, numAdjPolys, notInList, relationshipFound, sampleLineNum;
local array numeric adjPolyList[1];
local class POLYLINE currPlineOfNodes, tempPlineOfNodes;
local numeric temp1ClosestVertex, temp2ClosestVertex, dist1, dist2;
local class POINT2D sharedNodePt;
# While the there is a possibility for additional basins to be attached
numBasinsComputed = 1;
tempCount = 0;
statusC.BarInit(warnCountLostPolygons, 0);
while (numBasinsComputed > 0) {
tempCount += numBasinsComputed;
statusC.BarIncrement(numBasinsComputed, 0);
numBasinsComputed = 0;
# Find any basins without sample point ID's that should have one, and give them an ID
for i = 1 to numUserBasins {
numRecords = TableReadAttachment(UBasinPolyToPolyTBL, i, records, "polygon");
# If the polygon is not currently attached to any sample point
if (numRecords == 0) {
# Get all the lines for the polygon so the nodes can be determined
numNodes = GetVectorPolyLineList(UBasinVECT, lines1, i);
currPlineOfNodes.Clear();
for j = 1 to numNodes {
P = GetVectorLine(UBasinVECT, lines1[j]);
# Get the two end points and check if they already exist in the polyline of "nodes"
pt1 = P.GetVertex(0);
pt2 = P.GetVertex(P.GetNumPoints() - 1);
testVertex = currPlineOfNodes.FindClosestVertex(pt1);
testPt = currPlineOfNodes.GetVertex(testVertex);
# If it doesn't exist, add it
if (testPt != pt1) {
currPlineOfNodes.AppendVertex(pt1);
}
testVertex = currPlineOfNodes.FindClosestVertex(pt2);
testPt = currPlineOfNodes.GetVertex(testVertex);
# If it doesn't exist, add it
if (testPt != pt2) {
currPlineOfNodes.AppendVertex(pt2);
}
}
# Get all the adjacent polygons. Any connected only by a node (which is what the script is looking for)
# will be excluded in this list.
numAdjPolys = GetVectorPolyAdjacentPolyList(UBasinVECT, adjPolyList, i);
relationshipFound = 0;
# For all polgons
for j = 1 to numUserBasins {
if (relationshipFound == 0) {
# Make sure it doesn't check itself
if (j != i) {
# Make sure the polygon j is not in the list of adjacent polygons. Only looking for polygons that are
# connected by a single node, which are not classified as adjacent in the above function.
notInList = 1;
for k = 1 to numAdjPolys {
if ( j == adjPolyList[k]) {
notInList = 0;
}
}
if (notInList) {
# Get all the lines for the possibly "attached by node" polygon so the nodes can be determined
numNodes = GetVectorPolyLineList(UBasinVECT, lines2, j);
# For each line, get the endpoints and save them in a polyline
tempPlineOfNodes.Clear();
for k = 1 to numNodes {
P = GetVectorLine(UBasinVECT, lines2[k]);
# Get the two end points and check if they already exist in the polyline of "nodes"
pt1 = P.GetVertex(0);
pt2 = P.GetVertex(P.GetNumPoints() - 1);
testVertex = tempPlineOfNodes.FindClosestVertex(pt1);
testPt = tempPlineOfNodes.GetVertex(testVertex);
# If it doesn't exist, add it
if (testPt != pt1) {
tempPlineOfNodes.AppendVertex(pt1);
}
testVertex = tempPlineOfNodes.FindClosestVertex(pt2);
testPt = tempPlineOfNodes.GetVertex(testVertex);
# If it doesn't exist, add it
if (testPt != pt2) {
tempPlineOfNodes.AppendVertex(pt2);
}
}
# Compare all the "nodes" saved in each polyline.
for m = 0 to currPlineOfNodes.GetNumPoints() - 1 {
for n = 0 to tempPlineOfNodes.GetNumPoints() - 1 {
pt1 = currPlineOfNodes.GetVertex(m);
pt2 = tempPlineOfNodes.GetVertex(n);
if (pt1 == pt2) {
# Save a copy of the point for later test
sharedNodePt = pt1;
# Find the closest point on one of the flowpath lines
pt1 = ObjectToMap(UBasinVECT, pt1.x, pt1.y, UBSNgeoref);
lineNum = FindClosestLine(STDFlowPath, pt1.x, pt1.y, FLOWgeoref);
tempPt1 = MapToObject(FLOWgeoref, pt1.x, pt1.y, STDFlowPath);
ClosestPointOnLine(STDFlowPath, lineNum, tempPt1.x, tempPt1.y, tempValx, tempValy);
tempPt2.x = tempValx;
tempPt2.y = tempValy;
pt2 = ObjectToMap(STDFlowPath, tempPt2.x, tempPt2.y, FLOWgeoref);
# Find the corner of the raster cell that the point are located at
pt1 = MapToObject(DEMgeoref, pt1.x, pt1.y, SampleDEM);
pt2 = MapToObject(DEMgeoref, pt2.x, pt2.y, SampleDEM);
pt1.x = round(pt1.x);
pt1.y = round(pt1.y);
pt2.x = round(pt2.x);
pt2.y = round(pt2.y);
# If they flowpath line goes through the shared node, then the main polygon is either
# upstream or downstream of the current polygon.
if ((pt1.x == pt2.x) && (pt1.y == pt2.y)) {
# Make sure the flowpath that goes through the node is going through the polygons
# that are being checked. pt2 is still the shared node in. By setting a test point
# a half object cordinate in each diagnal direction and then snapping it to the
# flowpath line, it can then check if the flowpath runs through both of the polygons
# that are being tested (i and j). This will stop cases in which the polygons attached
# only by a shared node, do not have a flowpath that connects the two.
stillPassingTest = 1;
for t = 0 to 1 {
for u = 0 to 1 {
if (stillPassingTest == 1) {
# Set the test point in one of the 4 diagnal directions and then snap it to the flowpath
tempPt1.x = sharedNodePt.x - 0.5 + t;
tempPt1.y = sharedNodePt.y - 0.5 + u;
tempPt1 = ObjectToMap(UBasinVECT, tempPt1.x, tempPt1.y, UBSNgeoref);
tempPt1 = MapToObject(FLOWgeoref, tempPt1.x, tempPt1.y, STDFlowPath);
ClosestPointOnLine(STDFlowPath, lineNum, tempPt1.x, tempPt1.y, tempValx, tempValy);
tempPt1.x = tempValx;
tempPt1.y = tempValy;
tempPt1 = ObjectToMap(STDFlowPath, tempPt1.x, tempPt1.y, FLOWgeoref);
tempPt1 = MapToObject(UBSNgeoref, tempPt1.x, tempPt1.y, UBasinVECT);
# Reset to object coordinates in 0.5 divisions to compensate for minor rounding problems
tempPt1.x = round(tempPt1.x * 2) / 2;
tempPt1.y = round(tempPt1.y * 2) / 2;
# If the snapped test point is not the same as the shared point, find out which
# polygon it is in.
if (tempPt1 != sharedNodePt) {
tempPt1 = ObjectToMap(UBasinVECT, tempPt1.x, tempPt1.y, UBSNgeoref);
polyNum = FindClosestPoly(UBasinVECT, tempPt1.x, tempPt1.y, UBSNgeoref, 0);
# If it is in either of the two polygons that are being checked, it's ok
if ( (polyNum == i) or (polyNum == j) ) {
stillPassingTest = 1;
# If it's not in the two polygons, there is no relationship
} else {
stillPassingTest = 0;
}
# If the snapped test point is the same as the shared point, keep checking
} else {
stillPassingTest = 1;
}
}
}
}
if (stillPassingTest) {
# Check if the attached by node polygon has a sample ID number
numRecords = TableReadAttachment(UBasinPolyToPolyTBL, j, records, "polygon");
if (numRecords != 0) {
#if (sampleNum$.length != 0) {
sampleNum$ = TableReadFieldStr(UBasinPolyToPolyTBL, sampleIDFldName, records[1]);
# Find out what line the sample point is on
records[1] = TableKeyFieldLookup(ResultSampleTBL, sampleIDFldName, sampleNum$, "Equals");
TableGetRecordElementList(ResultSampleTBL, records[1], elements, "point");
TableReadAttachment(ResultPointToPointTBL, elements[1], records, "point");
sampleLineNum = TableReadFieldNum(ResultPointToPointTBL, onLineFldName, records[1]);
# If the both the shared node and the sample point are on the same flowpath line segment
if (sampleLineNum == lineNum) {
# Get the line and make sure the downstream node is the index of 0
P = GetVectorLine(STDFlowPath, lineNum);
# Get one of the endpoints
pt1 = P.GetVertex(0);
ObjectToMap(STDFlowPath, pt1.x, pt1.y, FLOWgeoref, Temp1MapX, Temp1MapY);
node1 = FindClosestNode(STDFlowPath, Temp1MapX, Temp1MapY, FLOWgeoref);
# Node 1 needs to be downstream. If it's not, turn the Polyline around
if (node1 != STDFlowPath.line[lineNum].(LineToLineTblName).(downNodeFldName)) {
P.Reverse();
}
# tempPt2 is the shared node point
# tempPt1 will be the sample point. Assign the values from the sample point.
records[1] = TableKeyFieldLookup(ResultSampleTBL, sampleIDFldName, sampleNum$, "Equals");
TableGetRecordElementList(ResultSampleTBL, records[1], elements, "point");
tempPt1.x = ResultVECT.point[elements[1]].Internal.x;
tempPt1.y = ResultVECT.point[elements[1]].Internal.y;
tempPt1 = ObjectToMap(ResultVECT, tempPt1.x, tempPt1.y, RSLTgeoref);
tempPt1 = MapToObject(FLOWgeoref, tempPt1.x, tempPt1.y, STDFlowPath);
temp1ClosestVertex = P.FindClosestVertex(tempPt1);
temp2ClosestVertex = P.FindClosestVertex(tempPt2);
# If they have the same vertex, then determine which is further downstream by tests
if (temp1ClosestVertex == temp2ClosestVertex) {
testVertex = temp1ClosestVertex - 1;
# If the tempVerticies are 0, they can only be upstream of it, so use it for the test
if (testVertex == -1) {
testVertex = 0;
}
# Find the distances to the downstream test point
testPt = P.GetVertex(testVertex);
Displacement3D(tempPt1.x, tempPt1.y, 0, testPt.x, testPt.y, 0, dist1, junk1, junk2);
Displacement3D(tempPt2.x, tempPt2.y, 0, testPt.x, testPt.y, 0, dist2, junk1, junk2);
# If the sample point is closer to the downstream vertex, then the shared node
# is upstream of the sample point.
if (dist1 < dist2) {
relationshipFound = j;
}
# If temp1's vertex less then temp2's, then it is downstream
} else if (temp1ClosestVertex < temp2ClosestVertex) {
relationshipFound = j;
# Otherwise it is not downstream
}
# If their not, go downstream and check if the line the sample point is on, is downstream
# of the line that the shared node is on.
} else {
while (lineNum != -1) {
TableReadAttachment(STDFlowPathTBL, lineNum, records, "line");
lineNum = TableReadFieldNum(STDFlowPathTBL, downLineFldName, records[1]);
# If the sample's line IS downstream of the shared node, the polygon should have the
# same sample ID as the polygon it shares a node with.
if (lineNum == sampleLineNum) {
relationshipFound = j;
lineNum = -1;
}
} # End While
} # End Vertex Comparison
} # End length != 0
} # End stillPassingTest
} # End if shared node is on flowpath
} # End if they share a node
} # End for n
} # End for m
} # End if not on list
} # End if not itself
} # End if relationship is not found
} # End for j loop to numBasins
# If a relationship was found, record it
if (relationshipFound) {
# Attach the polygon to the record that has the correct sample point ID
TableReadAttachment(UBasinPolyToPolyTBL, relationshipFound, records, "polygon");
TableWriteAttachment(UBasinPolyToPolyTBL, i, records, 1, "polygon");
# Add 1 to the number of polygons previously attached to the sample point
numPolys = TableReadFieldNum(UBasinPolyToPolyTBL, numPolysPerSampleFldName, records[1]) + 1;
TableWriteField(UBasinPolyToPolyTBL, records[1], numPolysPerSampleFldName, numPolys);
message = "Correction: The following polygon element " + NumToStr(i) + " does not have a sample point inside of\n" +
" it, but since it was connected to another polygon by only a single node and\n" +
" an intersecting flowpath line, it was attached to the sample point " + sampleNum$ + ".\n\n";
print(message);
fwritestring(summaryFile, message);
# Increment the counters
numBasinsComputed = numBasinsComputed + 1;
warnCountManAttachPolygons = warnCountManAttachPolygons + 1;
}
} # End if i doesn't have a sampleID
} # End for i loop
} # End While loop
message = "##########################################################################################\n" +
"Number of manually attached polygons found: " + NumToStr(warnCountManAttachPolygons) + "\n" +
"##########################################################################################\n\n";
print(message);
fwritestring(summaryFile, message);
#######################################################################################################
# Attach the user defined catchment basin polygons to the copied sediment sample table
statusC.Message = "Step 16 of 18: Attatching polygons to the sediment sample database";
#######################################################################################################
statusC.BarInit(numUserBasins, 0);
for i = 1 to numUserBasins {
statusC.BarIncrement(1, 0);
# If the polygon is attached to a sample in the PolyToPoly Table, transfer that attachment to the sample table
if (TableReadAttachment(UBasinPolyToPolyTBL, i, records, "polygon") > 0) {
sampleNum$ = TableReadFieldStr(UBasinPolyToPolyTBL, sampleIDFldName, records[1]);
# Make sure the basin has a sample point
if (sampleNum$.length != 0) {
records[1] = TableKeyFieldLookup(UBasinSampleTBL, sampleIDFldName, sampleNum$, "Equals");
TableWriteAttachment(UBasinSampleTBL, i, records, 1, "polygon");
}
}
}
#######################################################################################################
# Calculate the statistics for the combined upstream basins for each point
statusC.Message = "Step 17 of 18: Computing upstream area statistics for each catchment basin";
#######################################################################################################
TableAddFieldInteger(UBasinPolyToPolyTBL, totalSamplesIncludedFldName, 15); # Number of Samples included (self included)
TableAddFieldFloat(UBasinPolyToPolyTBL, totalSampleBasinAreaFldName, 15, 4); # Total area of all polygons upstream of sample point (self included)
TableAddFieldFloat(UBasinSampleTBL, totalSampleBasinAreaFldName, 15, 4); # Total area of all polygons upstream of sample point (self included)
TableAddFieldInteger(UBasinPolyToPolyTBL, totalSampleBasinsFldName, 12); # Total number of sample basins for sample point (self included)
TableAddFieldInteger(UBasinPolyToPolyTBL, numUpBasinsFldName, 12); # Number of basins upstream (self excluded)
local numeric numUpSamples, totalUpSamples, tempAreaValue, numElements;
local numeric totalInclArea, currTotalArea, upstreamSamplesTotalArea;
if (DEBUG) {
statusC.Message = "Step 17 of 18: Computing upstream area statistics for each catchment basin (FIRST LOOP)";
fwritestring(summaryFile, "Step 17 of 18: Computing upstream area statistics for each catchment basin (FIRST LOOP)" + "\n");
# PopupMessage("Start debugging here");
}
# Determine the area of the basin for each sample point. This will be the sum of all basins with
# the same sample point ID
statusC.BarInit(numUserBasins, 0);
for i = 1 to numUserBasins {
statusC.BarIncrement(1, 0);
tempnum = TableReadAttachment(UBasinPolyToPolyTBL, i, records, "polygon");
sampleNum$ = TableReadFieldStr(UBasinPolyToPolyTBL, sampleIDFldName, records[1]);
if (DEBUG) {
fwritestring(summaryFile, " i:" + NumToStr(i) + " numS: " + sampleNum$ + " tempN:" + NumToStr(tempnum) + " rec[1]:" + NumToStr(records[1]) + "\n");
}
# Make sure the basin has a sample point
if (sampleNum$.length != 0) {
# Check if the basin has it's sample point area computed yet
tempAreaValue = TableReadFieldNum(UBasinPolyToPolyTBL, areaSamplePolysFldName, records[1]);
if (DEBUG) {
fwritestring(summaryFile, " tempA:" + NumToStr(tempAreaValue) + "\n");
}
# If it is invalid, it needs to compute it
if (tempAreaValue == -1) {
# Find common record for the polygons attached to the current sample point
records[1] = TableKeyFieldLookup(UBasinPolyToPolyTBL, sampleIDFldName, sampleNum$, "Equals");
if (DEBUG) {
fwritestring(summaryFile, " rec[1]:" + NumToStr(records[1]) + "\n");
}
tempAreaValue = 0;
# For all the polygons attached to that sample point, sum up their areas in square km
numElements = TableGetRecordElementList(UBasinPolyToPolyTBL, records[1], elements, "polygon");
for j = 1 to numElements {
tempAreaValue = tempAreaValue + (UBasinVECT.poly[elements[j]].POLYSTATS.Area / 1000000);
}
if (DEBUG) {
fwritestring(summaryFile, " tempA:" + NumToStr(tempAreaValue) + " numE:" + NumToStr(numElements) + "\n");
}
# Record the value back to the record for the polygons
TableWriteField(UBasinPolyToPolyTBL, records[1], areaSamplePolysFldName, tempAreaValue);
}
}
}
numBasinsComputed = 0;
if (DEBUG) {
statusC.Message = "Step 17 of 18: Computing upstream area statistics for each catchment basin (SECOND LOOP)";
fwritestring(summaryFile, "Step 17 of 18: Computing upstream area statistics for each catchment basin (SECOND LOOP)" + "\n");
}
# Copy the number of upstream samples for each sample to the polygon table, same as number of basins
statusC.BarInit(numSamplePoints, 0);
for i = 1 to numSamplePoints {
statusC.BarIncrement(1, 0);
# Get the sample name and number of upstream samples for the current sample point
TableReadAttachment(ResultSampleTBL, i, records, "point");
sampleNum$ = TableReadFieldStr(ResultSampleTBL, sampleIDFldName, records[1]);
TableReadAttachment(ResultPointToPointTBL, i, records, "point");
totalUpSamples = TableReadFieldNum(ResultPointToPointTBL, totalUpSamplesFldName, records[1]);
numUpSamples = TableReadFieldNum(ResultPointToPointTBL, numUpSamplesFldName, records[1]);
# Find the record and copy some of the fields from the points table to the polygon table
records[1] = TableKeyFieldLookup(UBasinPolyToPolyTBL, sampleIDFldName, sampleNum$, "Equals");
TableWriteField(UBasinPolyToPolyTBL, records[1], totalSamplesIncludedFldName, totalUpSamples + 1);
TableWriteField(UBasinPolyToPolyTBL, records[1], totalSampleBasinsFldName, totalUpSamples);
TableWriteField(UBasinPolyToPolyTBL, records[1], numUpBasinsFldName, numUpSamples);
# Start the process of computing the area for the overall sample basin by computing the basins
# with no upstream sample points. Use the area value of the sample basin for the "total"
# upstream basin value. Check if the sample basin has no upstream sample points.
if (numUpSamples == 0) {
# Assign the total area, the area of the sample basin
tempAreaValue = TableReadFieldNum(UBasinPolyToPolyTBL, areaSamplePolysFldName, records[1]);
TableWriteField(UBasinPolyToPolyTBL, records[1], totalSampleBasinAreaFldName, tempAreaValue);
numBasinsComputed = numBasinsComputed + 1;
# Otherwise, set the total number of upstream basins and the area to invalid
} else {
TableWriteField(UBasinPolyToPolyTBL, records[1], totalSampleBasinAreaFldName, -1);
}
}
tempCount = 0;
# Atleast 1 basin has to be the most upstream, so it should always enter the while loop.
statusC.BarInit(numUserBasins, 0);
while (numBasinsComputed > 0) {
if (DEBUG) {
fwritestring(summaryFile, " numB: " + NumToStr(numBasinsComputed) + " numS: " + NumToStr(numSamplePoints) + "\n");
}
# Reset the counter so the while loop will only stop when no new basin values are computed
tempCount += numBasinsComputed;
statusC.BarIncrement(numBasinsComputed, 0);
numBasinsComputed = 0;
for i = 1 to numSamplePoints {
# Get the polygon element number for the sample point
TableReadAttachment(ResultSampleTBL, i, records, "point");
sampleNum$ = TableReadFieldStr(ResultSampleTBL, sampleIDFldName, records[1]);
records[1] = TableKeyFieldLookup(UBasinPolyToPolyTBL, sampleIDFldName, sampleNum$, "Equals");
TableGetRecordElementList(UBasinPolyToPolyTBL, records[1], elements, "polygon");
# Get the current total included area for that basin
TableReadAttachment(UBasinPolyToPolyTBL, elements[1], records, "polygon");
totalInclArea = TableReadFieldNum(UBasinPolyToPolyTBL, totalSampleBasinAreaFldName, records[1]);
# If the total area is still invalid, we need to try and compute it
if (totalInclArea == -1) {
# Add the area for the current polygon into the total
currTotalArea = UBasinVECT.Poly[elements[1]].POLYSTATS.Area / 1000000;
# Get the number of upstream samples
TableReadAttachment(ResultPointToPointTBL, i, records, "point");
numUpstream = TableReadFieldNum(ResultPointToPointTBL, numUpSamplesFldName, records[1]);
# Initialize the list and then get the list of upstream sample numbers
UpSample.Clear();
UpSample.AddToFront("");
for j = 1 to numUpstream {
UpSample.AddToEnd( TableReadFieldStr(ResultPointToPointTBL, upSampleFldName + NumToStr(j), records[1]) );
}
if (DEBUG) {
fwritestring(summaryFile, " sNum: " + sampleNum$ + " numS: " + NumToStr(numSamplePoints) + " i:" + NumToStr(i) + " numU: " + NumToStr(numUpstream) + "\n");
}
for j = 1 to numUpstream {
# If the current total is not invalid, keep checking the upstream samples
if ( currTotalArea != -1) {
# Get the polygon element that the Sample number refers to and check if it's total area is valid
records[1] = TableKeyFieldLookup(UBasinPolyToPolyTBL, sampleIDFldName, UpSample[j], "Equals");
TableGetRecordElementList(UBasinPolyToPolyTBL, records[1], elements, "polygon");
TableReadAttachment(UBasinPolyToPolyTBL, elements[1], records, "polygon");
upstreamSamplesTotalArea = TableReadFieldNum(UBasinPolyToPolyTBL, totalSampleBasinAreaFldName, records[1]);
# If the upstream basin of sample point i has a total basin area value, add it to the current total
if (upstreamSamplesTotalArea != -1) {
currTotalArea = currTotalArea + upstreamSamplesTotalArea;
# Otherwise, the it hasn't been calculated, so set the current total to invalid to stop checking
} else {
currTotalArea = -1;
}
}
} # End for j loop
# If the total included area that was computed is valid, record it to the polygon's table
if ( currTotalArea != -1 ) {
# Get the polygon element number for the sample point
TableReadAttachment(ResultSampleTBL, i, records, "point");
sampleNum$ = TableReadFieldStr(ResultSampleTBL, sampleIDFldName, records[1]);
numRecords = TableKeyFieldLookupList(UBasinPolyToPolyTBL, sampleIDFldName, sampleNum$, records, "Equals");
# Write the new total area to the polygon table
for j = 1 to numRecords {
TableWriteField(UBasinPolyToPolyTBL, records[j], totalSampleBasinAreaFldName, currTotalArea);
}
# Update the while counter so it will run again
numBasinsComputed = numBasinsComputed + 1;
if (DEBUG) {
fwritestring(summaryFile, " currT: " + NumToStr(currTotalArea) + " sNum:" + sampleNum$ + " numR:" + NumToStr(numRecords) + "\n");
}
}
} # End if not invalid
} # End for i loop
} # End while loop
if (DEBUG) {
statusC.Message = "Step 17 of 18: Computing upstream area statistics for each catchment basin (THIRD LOOP)";
fwritestring(summaryFile, "Step 17 of 18: Computing upstream area statistics for each catchment basin (THIRD LOOP)" + "\n");
}
# Copy the Upstream Catchment Basin values to the point and polygon sediment sample databases
statusC.BarInit(numUserBasins, 0);
for i = 1 to numUserBasins {
statusC.BarIncrement(1, 0);
TableReadAttachment(UBasinPolyToPolyTBL, i, records, "polygon");
sampleNum$ = TableReadFieldStr(UBasinPolyToPolyTBL, sampleIDFldName, records[1]);
# Make sure the basin has a sample point
if (sampleNum$.length != 0) {
tempAreaValue = TableReadFieldNum(UBasinPolyToPolyTBL, totalSampleBasinAreaFldName, records[1]);
if (DEBUG) {
fwritestring(summaryFile, "i: " + NumToStr(i) + " tempA:" + NumToStr(tempAreaValue) + "\n");
}
# Find the record with the same SampleID and write the total area value to it.
# The same record may be written two twice if multiple polygons are attached to it, but the value should be the same.
records[1] = TableKeyFieldLookup(ResultSampleTBL, sampleIDFldName, sampleNum$, "Equals");
TableWriteField(ResultSampleTBL, records[1], totalSampleBasinAreaFldName, tempAreaValue);
records[1] = TableKeyFieldLookup(UBasinSampleTBL, sampleIDFldName, sampleNum$, "Equals");
TableWriteField(UBasinSampleTBL, records[1], totalSampleBasinAreaFldName, tempAreaValue);
}
}
#######################################################################################################
# Create a user defined upstream flowpath vector (clipped STDFLOWPATH inside USDBASIN polygons)
statusC.Message = "Step 18 of 18: Computing the user defined upstream flowpath vector";
statusC.BarClear(0);
#######################################################################################################
if ( mainDLG.GetCtrlByID("id_COMP_USR_UP_FLOWPATH").GetValueNum()) {
class VECTOR USDUpFlowVECT;
CreateVector(USDUpFlowVECT, OutputProjectFile, "USD_UPFLOWPATH", "Upstream flowpath from sample points");
USDUpFlowVECT = VectorExtract(UBasinVECT, STDFlowPath, "InsideClip", "None", "1", "1", "1");
CloseVector(USDUpFlowVECT);
}
#######################################################################################################
# End watershed main process
#######################################################################################################
CloseRaster(SampleDEM);
CloseVector(SampleVECT);
CloseVector(WarpVECT);
CloseVector(ResultVECT);
CloseVector(STDFlowPath);
CloseVector(UBasinVECT);
}
############
