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