######################################################################################################
#
# SampleCatchments.sml (Stand Alone)
#
######################################################################################################
#
# Author:
# Dave Breitwisch, MicroImages, Inc.
#
# Started:
# 19 September 2005
#
# Updated:
# 9 December 2011 - Fixed problem evaluating unattached basins in series
# 21 November 2011 - Increased size of array used while evaluating unattached basins
# 11 November 2011 - Fix checking flowdirs raster to fix unattached basins
#
# 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 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;
class 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 TRANS2D_MAPGEN 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$ = '<?xml version="1.0"?>
<root>
<dialog id = "MAIN_DIALOG" title="Multi-Sample Watershed Process and Analysis" Buttons = "">
<groupbox Name = "Select Input Data" ExtraBorder = "4">
<pane Orientation = "Horizontal">
<pane Orientation = "Vertical">
<pushbutton id = "id_INPUT_DEM_BUTTON" Name = "Input DEM... " VertResize = "Fixed" OnPressed = "OnGetInputDEMPressed();" Enabled="true"/>
<pushbutton id = "id_INPUT_VECTOR_BUTTON" Name = "Input Vector..." VertResize = "Fixed" OnPressed = "OnGetInputVECTORPressed();" Enabled="false"/>
</pane>
<pane Orientation = "Vertical">
<edittext id = "id_INPUT_DEM" Width = "70" ReadOnly = "true" Justify = "left"/>
<groupbox ExtraBorder = "2">
<edittext id = "id_INPUT_VECTOR" Width = "70" ReadOnly = "true" Justify = "left"/>
<pane Orientation = "Horizontal">
<label Width = "12" TextAlign = "right">Sample Table Name</label>
<edittext id = "id_TABLE_NAME" HorizResize = "Fixed" Width = "20" ReadOnly = "true"/>
<label Width = "12" TextAlign = "right">Unique ID Field</label>
<edittext id = "id_FIELD_NAME" HorizResize = "Fixed" Width = "20" ReadOnly = "true"/>
</pane>
</groupbox>
</pane>
</pane>
</groupbox>
<label/>
<groupbox Name = "Watershed Options">
<book>
<page Name = "General">
<pane Orientation = "Horizontal">
<pane Orientation = "Vertical" VertResize = "Fixed">
<togglebutton id = "id_FILL_ALL_DEP" Name = "Fill All Depressions" Selected = "true" Enabled = "true"/>
<togglebutton Name = "Compute Databases" Selected = "false" Enabled = "false"/>
</pane>
<pane Orientation = "Vertical" VertResize = "Fixed">
<togglebutton Name = "Compute Standard Flowpath" Selected = "true" Enabled = "false"/>
<togglebutton id = "id_COMP_STD_BASINS" Name = "Compute Standard Basins" Selected = "false" Enabled = "true"/>
<togglebutton id = "id_COMP_STD_RIDGES" Name = "Compute Standard Ridges" Selected = "false" Enabled = "true"/>
</pane>
<pane Orientation = "Vertical" VertResize = "Fixed">
<togglebutton Name = "Compute User-Defined Basins" Selected = "true" Enabled = "false"/>
<togglebutton id = "id_COMP_USR_FLOWPATH" Name = "Compute User-Defined Flowpath" Selected = "false" Enabled = "true"/>
<togglebutton id = "id_COMP_USR_UP_FLOWPATH" Name = "Compute User-Defined Upstream Flowpath" Selected = "true" Enabled = "true"/>
</pane>
</pane>
</page>
<page Name = "Mask">
<togglebutton Name = "Use Watershed Mask" id = "id_MASK_TOGGLE" Enabled = "true" Selected = "false" OnChanged = "OnUseMaskToggleChanged();"/>
<groupbox ExtraBorder = "2" VertResize = "Fixed">
<pane Orientation = "Horizontal" VertResize = "Fixed">
<pane Orientation = "Vertical" HorizResize = "Fixed">
<pushbutton id = "id_MASK_BUTTON" Name = "Mask..." HorizResize = "Fixed" Enabled = "false" OnPressed = "OnMaskButtonPressed();"/>
</pane>
<pane Orientation = "Vertical" HorizResize = "Expand">
<edittext id = "id_INPUT_MASK" HorizResize = "Expand" Enabled = "false" ReadOnly = "true"/>
</pane>
</pane>
</groupbox>
</page>
<page Name = "Flow Path and Basin">
<pane HorizResize = "Fixed">
<pane Orientation = "Horizontal">
<label HorizResize = "Fixed">Threshold</label>
<pushbutton Name = "Restore Defaults" HorizResize = "Fixed" OnPressed = "OnRestoreDefaultsPressed();"/>
<togglebutton id = "id_SEP_VALLEY_POLYS" Name = "Separate Valley Polygons"/>
</pane>
<groupbox ExtraBorder = "2">
<pane Orientation = "Horizontal">
<pane Orientation = "Vertical" HorizResize = "Fixed">
<label>Outlet</label>
<label>Branch</label>
</pane>
<pane Orientation = "Vertical" HorizResize = "Expand">
<editnumber id = "id_OUTLET" Default = "256" Precision = "0"/>
<editnumber id = "id_BRANCH" Default = "128" Precision = "0"/>
</pane>
<pane Orientation = "Vertical" HorizResize = "Fixed">
<label>Inlet</label>
<label>Basin</label>
</pane>
<pane Orientation = "Vertical" HorizResize = "Expand">
<editnumber id = "id_INLET" Default = "32" Precision = "0"/>
<editnumber id = "id_BASIN" Default = "378" Precision = "0"/>
</pane>
</pane>
</groupbox>
</pane>
</page>
</book>
</groupbox>
<label/>
<groupbox Name = "Sample Location Adjustment Options" ExtraBorder = "4">
<pane Orientation = "Horizontal">
<groupbox ExtraBorder = "2">
<pane Orientation = "Vertical">
<togglebutton id = "id_SNAP_LINE_TOGGLE" Name = "Snap to closest flowpath" Selected = "true" Enabled = "false" OnChanged = "OnSnapLineToggleChanged();"/>
<pane Orientation = "Horizontal">
<label id = "id_SNAP_LINE_LABEL" TextAlign = "Right">Maximum search distance (meters)</label>
<editnumber id = "id_SNAP_LINE_MAX_DIST" Width = "10" Precision = "2" Default = "200"/>
</pane>
</pane>
</groupbox>
<groupbox ExtraBorder = "2">
<pane Orientation = "Vertical">
<togglebutton id = "id_SNAP_NODE_TOGGLE" Name = "Snap to downstream node" Selected = "false" OnChanged = "OnSnapNodeToggleChanged();"/>
<pane Orientation = "Horizontal">
<label id = "id_SNAP_NODE_LABEL" TextAlign = "Right" Enabled = "false">Maximum search distance (meters)</label>
<editnumber id = "id_SNAP_NODE_MAX_DIST" Width = "10" Precision = "2" Enabled = "false" Default = "90"/>
</pane>
</pane>
</groupbox>
</pane>
</groupbox>
<label/>
<groupbox Name = "Select Output Locations" ExtraBorder = "4">
<pane Orientation = "Horizontal">
<pane Orientation = "Vertical">
<pushbutton id = "id_OUTPUT_PROJECT_FILE_BUTTON" Name = "Output Project File..." OnPressed = "OnGetOutputProjectFilePressed();" Enabled="false"/>
<pushbutton id = "id_OUTPUT_SUMMARY_FILE_BUTTON" Name = "Output Summary File..." OnPressed = "OnGetOutputSummaryFilePressed();" Enabled="false"/>
</pane>
<pane Orientation = "Vertical">
<edittext id = "id_OUTPUT_PROJECT_FILE" Width = "70" ReadOnly = "true" Justify = "left"/>
<edittext id = "id_OUTPUT_SUMMARY_FILE" Width = "70" ReadOnly = "true" Justify = "left"/>
</pane>
</pane>
</groupbox>
<pane Orientation = "Horizontal" HorizResize = "Expand">
<pane Orientation = "Vertical"/>
<pushbutton id = "id_RUN_BUTTON" Name = "Run..." HorizResize = "Fixed" Enabled = "false" OnPressed = "OnRunPressed();"/>
<pushbutton id = "id_CANCEL_BUTTON" Name = "Cancel" HorizResize = "Fixed" Enabled = "1" OnPressed = "OnCancelPressed();"/>
</pane>
</dialog>
</root>';
######################################################################################################
######################################################################################################
######################################################################################################
# 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(index, 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" +
"Warning: 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 numeric 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 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 closest 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 sample 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 related 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 its 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);
#######################################################################################################
# Iterate through basin polygons that have no attached sample record; trace flow directions downstream from a
# cell inside the basin to find if the basin connects to an attributed downstream basin. If so, assign
# the associated sample record to the current polygon.
# This solves the case in which a catchment basin is "pinched" due to a diagonal
# flowpath line, creating separate basin polygons connected only by a node.
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
WatershedGetObject(watershed, "RasterFlowDirection", ObjectFilename, ObjectName);
class RVC_RASTER FlowDirs;
OpenRaster(FlowDirs, ObjectFilename, ObjectName);
numeric nullval = NullValue(FlowDirs);
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) {
array numeric UnattachedPolys[warnCountLostPolygons];
numeric index = 1;
for index = 1 to warnCountLostPolygons {
UnattachedPolys[index] = -1;
}
index = 1;
UnattachedPolys[index] = i;
relationshipFound = 0;
class POLYLINE poly = GetVectorPoly(UBasinVECT, i);
class POINT2D vertex = poly.GetVertex(1);
vertex.x += .5;
vertex.y -= .5;
if (!poly.IsPointInside(vertex)) {
vertex = poly.GetVertex(1);
vertex.x -= .5;
vertex.y -= .5;
}
if (!poly.IsPointInside(vertex)) {
vertex = poly.GetVertex(1);
vertex.x -= .5;
vertex.y += .5;
}
if (!poly.IsPointInside(vertex)) {
vertex = poly.GetVertex(1);
vertex.x += .5;
vertex.y += .5;
}
while (true) {
class POINT2D FlowCell = ObjectToMap(UBasinVECT, vertex.x, vertex.y, UBSNgeoref);
FlowCell = MapToObject(UBSNgeoref, FlowCell.x, FlowCell.y, FlowDirs);
numeric FlowValue = FlowDirs[floor(FlowCell.y), floor(FlowCell.x)];
if (FlowValue == nullval) break;
if (IsNull(FlowValue)) break;
switch (FlowValue) {
case 1: # upper right
vertex.x += 1;
vertex.y += 1;
break;
case 2: # right
vertex.x += 1;
break;
case 3: # lower right
vertex.x += 1;
vertex.y -= 1;
break;
case 4: # lower
vertex.y -= 1;
break;
case 5: # lower left
vertex.x -= 1;
vertex.y -= 1;
break;
case 6: # left
vertex.x -= 1;
break;
case 7: # upper left
vertex.x -= 1;
vertex.y += 1;
break;
case 8: # upper
vertex.y += 1;
break;
}
numeric AdjPoly = FindClosestPoly(UBasinVECT, vertex.x, vertex.y);
class POLYLINE AdjPolyline = GetVectorPoly(UBasinVECT, AdjPoly);
if (AdjPoly == 0 || !AdjPolyline.IsPointInside(vertex)) {
break;
}
if (AdjPoly != UnattachedPolys[index]) {
numRecords = TableReadAttachment(UBasinPolyToPolyTBL, AdjPoly, records, "polygon");
if (numRecords > 0) {
relationshipFound = AdjPoly;
break;
}
else {
UnattachedPolys[++index] = AdjPoly;
}
}
}
# 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");
for index = 1 to 100 {
if (UnattachedPolys[index] == -1) break;
TableWriteAttachment(UBasinPolyToPolyTBL, UnattachedPolys[index], 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(UnattachedPolys[index]) + " does not have a sample point inside of\n" +
" it, but since it was connected to an attributed downstream polygon, it was attached to that sample point.\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);
}
######################################################################################################
# Dialog Related Functions
######################################################################################################
######################################################################################################
# Called when "Run..." is pressed.
proc OnRunPressed( ) {
mainDLG.GetCtrlByID("id_RUN_BUTTON").SetEnabled(0);
if (sampleTblName == "" or sampleIDFldName == "") {
message = "The unique sample ID table and field values are not selected.\n";
message += "Please reselect the input vector and select the unique table and field";
PopupMessage(message);
mainDLG.GetCtrlByID("id_RUN_BUTTON").SetEnabled(1);
print(message);
fwritestring(summaryFile, message);
return;
}
# Set up the summary file
summaryFile = fopen(OutputSummaryFile, "w");
# Call Starting Procedure Here!!!!!
statusD.Create();
statusC = statusD.CreateContext();
string test = SampleIDThatIsNotUnique();
if ( test != "" ) {
InputVECTORfilename = InputVECTORobjectName = sampleTblName = sampleIDFldName = "";
mainDLG.GetCtrlByID("id_INPUT_VECTOR").SetValueStr("");
mainDLG.GetCtrlByID("id_TABLE_NAME").SetValueStr("");
mainDLG.GetCtrlByID("id_FIELD_NAME").SetValueStr("");
statusD.destroy();
message = "The following sample ID in the selected table is not unique: " + test + "\n";
message += "Please resolve the situation in the table before running the script again, \n";
message += "or select a new table and field by reselecting the input vector.";
PopupMessage(message);
mainDLG.GetCtrlByID("id_RUN_BUTTON").SetEnabled(1);
return;
}
# Set up a timer to record completion time
local class TIMER time;
time.Running = 0;
time.Value = 0;
time.Running = 1;
message = "Starting Watershed Process...\n\n";
print(message);
fwritestring(summaryFile, message);
StartWatershedProcess();
statusD.Destroy();
time.Running = 0;
message = sprintf("Time to complete: %2d minutes %0.2d seconds\n", floor(time.Value / 60), floor(time.Value % 60));
PopupMessage(message);
print(message);
fwritestring(summaryFile, message);
mainDLG.Close(1);
}
######################################################################################################
# Called when "Canel" is pressed. Closes dialog and exits script.
proc OnCancelPressed( ) {
mainDLG.Close(1);
}
######################################################################################################
# Called when "Input DEM..." is pressed. Opens dialog to select raster.
proc OnGetInputDEMPressed( ) {
local class RASTER DEM;
# Get the raster and set global variables
GetInputRaster(DEM);
class RVC_OBJITEM DEMobjitem = DEM.GetObjItem();
InputDEMfilename = DEMobjitem.GetFilePath();
InputDEMobjectName = DEMobjitem.GetObjectPath();
CloseRaster(DEM);
class FILEPATH demFilepath(InputDEMfilename);
# Update the dialog
local class GUI_CTRL_EDIT_STRING text = mainDLG.GetCtrlByID("id_INPUT_DEM");
text.SetValue(InputDEMfilename + " / " + InputDEMobjectName, 0);
mainDLG.GetCtrlByID("id_INPUT_VECTOR_BUTTON").SetEnabled(1);
}
######################################################################################################
# Called when "Input Vector" is pressed. Opens dialog to select vector.
proc OnGetInputVECTORPressed( ) {
local class VECTOR Vect;
# Get the Vector and set global variables
GetInputVector(Vect);
InputVECTORfilename = Vect.$Info.Filename;
InputVECTORobjectName = Vect.$Info.Name;
# Update the dialog
local class GUI_CTRL_EDIT_STRING text = mainDLG.GetCtrlByID("id_INPUT_VECTOR");
text.SetValue(InputVECTORfilename + " / " + InputVECTORobjectName, 0);
mainDLG.GetCtrlByID("id_OUTPUT_PROJECT_FILE_BUTTON").SetEnabled(1);
PopupSelectTableField(OpenVectorPointDatabase(Vect), sampleTblName, sampleIDFldName);
mainDLG.GetCtrlByID("id_TABLE_NAME").SetValueStr(sampleTblName);
mainDLG.GetCtrlByID("id_FIELD_NAME").SetValueStr(sampleIDFldName);
CloseVector(Vect);
}
######################################################################################################
# Called when the "Use Mask" toggle is pressed
proc OnUseMaskToggleChanged() {
numeric state = mainDLG.GetCtrlByID("id_MASK_TOGGLE").GetValueNum();
mainDLG.GetCtrlByID("id_MASK_BUTTON").SetEnabled(state);
mainDLG.GetCtrlByID("id_INPUT_MASK").SetEnabled(state);
}
######################################################################################################
# Called when "Mask..." is pressed.
proc OnMaskButtonPressed() {
local class RASTER Mask;
# Get the raster and set global variables
GetInputRaster(Mask);
InputMaskfilename = Mask.$Info.Filename;
InputMaskobjectName = Mask.$Info.Name;
CloseRaster(Mask);
# Update the dialog
local class GUI_CTRL_EDIT_STRING text = mainDLG.GetCtrlByID("id_INPUT_MASK");
text.SetValue(InputMaskfilename + " / " + InputMaskobjectName, 0);
}
######################################################################################################
# Called when "Restore Defaults" is pressed.
proc OnRestoreDefaultsPressed() {
mainDLG.GetCtrlByID("id_OUTLET").SetValueNum(256);
mainDLG.GetCtrlByID("id_BRANCH").SetValueNum(128);
mainDLG.GetCtrlByID("id_INLET").SetValueNum(32);
mainDLG.GetCtrlByID("id_BASIN").SetValueNum(378);
}
######################################################################################################
# Called when the "Snap to closest flowpath" toggle button is pressed
proc OnSnapLineToggleChanged() {
local class GUI_CTRL_TOGGLEBUTTON toggle = mainDLG.GetCtrlByID("id_SNAP_LINE_TOGGLE");
mainDLG.GetCtrlByID("id_SNAP_LINE_LABEL").SetEnabled(toggle.GetValue());
mainDLG.GetCtrlByID("id_SNAP_LINE_MAX_DIST").SetEnabled(toggle.GetValue());
}
######################################################################################################
# Called when "Snap to downstream node" toggle button is pressed
proc OnSnapNodeToggleChanged() {
local class GUI_CTRL_TOGGLEBUTTON toggle = mainDLG.GetCtrlByID("id_SNAP_NODE_TOGGLE");
mainDLG.GetCtrlByID("id_SNAP_NODE_LABEL").SetEnabled(toggle.GetValue());
mainDLG.GetCtrlByID("id_SNAP_NODE_MAX_DIST").SetEnabled(toggle.GetValue());
}
######################################################################################################
# Called when "Output Project File..." is pressed.
proc OnGetOutputProjectFilePressed( ) {
local class GUI_CTRL_EDIT_STRING text = mainDLG.GetCtrlByID("id_OUTPUT_PROJECT_FILE");
local class STRING tempStartDir = _context.ScriptDir;
if (InputDEMfilename != "") {
demFilepath.StripToFolder();
tempStartDir = demFilepath;
}
OutputProjectFile = GetOutputFileName(tempStartDir, "Select project file to save files to:", "rvc");
class FILEPATH outputFilepath(OutputProjectFile);
# Check if the user pressed "Cancel" in the dialog
if (OutputProjectFile == tempStartDir) {
text.ClearValue(0);
text.SetEnabled(0);
return;
}
# Delete the file if it exists. The user would have received a "Overwrite this file?" warning.
if (fexists(OutputProjectFile)) {
DeleteFile(OutputProjectFile);
}
CreateProjectFile(OutputProjectFile, "Created by SML Watershed Script to save resulting files");
text.SetValue(OutputProjectFile, 0);
mainDLG.GetCtrlByID("id_OUTPUT_SUMMARY_FILE_BUTTON").SetEnabled(1);
}
######################################################################################################
# Called when the "Get Summary File..." button is pressed.
proc OnGetOutputSummaryFilePressed() {
local class GUI_CTRL_EDIT_STRING text = mainDLG.GetCtrlByID("id_OUTPUT_SUMMARY_FILE");
outputFilepath.StripToFolder();
local class STRING defaultName$ = outputFilepath + "/Summary";
OutputSummaryFile = GetOutputFileName(defaultName$, "Select output text file", "txt");
text.SetValue(OutputSummaryFile, 0);
mainDLG.GetCtrlByID("id_RUN_BUTTON").SetEnabled(1);
}
######################################################################################################
# Procedure called to set up the dialog using the XML code specified in this script. Dialog is made
# using a Modal function call. When the dialog closes, the script will end.
proc OnCreateDialog( ) {
class XMLDOC mainDOC;
class XMLNODE mainNODE;
errXML = mainDOC.Parse(xmlCODE$);
if ( errXML < 0 ) {
PopupError( errXML ); # pop up an error dialog
Exit( );
}
mainNODE = mainDOC.GetElementByID( "MAIN_DIALOG" );
if ( mainNODE == 0 ) {
PopupMessage( "Could not find main dialog node in XML document" );
Exit();
}
mainDLG.SetXMLNode(mainNODE);
errDLG = mainDLG.DoModal();
if ( errDLG == -1 ) {
Exit();
}
}
######################################################################################################
Download