Accuracy evaluation of a 3-D spatial modeling approach to model linear objects and predict their lengths

target point 'and determine the centerpoints for these surrounding cells Dim I1 As Integer Dim I2 As Integer Dim D1 As Double Dim D2 As Double 'pIDArray and pRIDObj are used for the target point while the others are used for the surrounding cell centers. Dim pIDArray As IArray Dim pRIDObj As IRasterIdentifyObj Dim pIDArray1 As IArray Dim pRIDObj1 As IRasterIdentifyObj Dim pIDArray2 As IArray Dim pRIDObj2 As IRasterIdentifyObj Dim pIDArray3 As IArray Dim pRIDObj3 As IRasterIdentifyObj Dim pIDArray4 As IArray Dim pRIDObj4 As IRasterIdentifyObj Dim i As Long Dim j As Long Dim pPoint As IPoint Dim pFeature As IFeature 'pNewPoint is the target point, pNewPoint1, 2, 3, 4 are the center points of the surrounding cells. Dim pNewPoint As IPoint Set pNewPoint = New Point Dim pNewPoint1 As IPoint Set pNewPoint1 = New Point Dim pNewPoint2 As IPoint Set pNewPoint2 = New Point Dim pNewPoint3 As IPoint Set pNewPoint3 = New Point Dim pNewPoint4 As IPoint Set pNewPoint4 = New Point 'Loop through each line in the feature class and obtain points 'including start point and end point and points with cellsize interval Dim NumOfRow As Integer NumOfRow = pInFeatureClass.FeatureCount (Nothing) Dim pLine As IPolyline Set pLine = New Polyline Dim length2d As Double For i = 0 To NumOfRow - 1 'Get line Set pFeature = pInFeatureClass.GetFeature (i) Set pLine = pFeature.Shape length2d = pLine.Length 'jj and jjj will be used to determine how many points will be obtained Dim jj As Double Dim jjj As Integer Dim interval As Integer interval = cellwidth / 2 430 jj = length2d / interval jjj = length2d / interval If (jj - jjj) > 0 Then jjj = jjj + 1 Else End If For j = 0 To jjj Dim myrow As IRow Set myrow = pIntable.CreateRow myrow.value (5) = pFeature.value (2) myrow.value (1) = pLine.Length If j = 0 Then 'Get the start point Set pPoint = pLine.FromPoint myrow.value (2) = 0 myrow.value (3) = pLine.FromPoint.X myrow.value (4) = pLine.FromPoint.Y ElseIf j >= jjj Then 'Get the end point Set pPoint = pLine.ToPoint myrow.value (2) = pLine.Length myrow.value (3) = pLine.ToPoint.X myrow.value (4) = pLine.ToPoint.Y Else 'Get intermediate points based on the given interval pLine.QueryPoint esriNoExtension, j * interval, False, pPoint myrow.value (2) = j * interval myrow.value (3) = pPoint.X myrow.value (4) = pPoint.Y End If pNewPoint.X = pPoint.X pNewPoint.Y = pPoint.Y 'Get RasterIdentifyObject on that point Set pIDArray = pIdentify.Identify (pNewPoint) I1 = (pNewPoint.X - XMIN) / cellwidth I2 = (YMAX - pNewPoint.Y) / cellheight 'Identify the center points of the surrounding cells pNewPoint1.X = XMIN + (I1 - 0.5) * cellwidth pNewPoint1.Y = YMAX - (I2 - 0.5) * cellheight pNewPoint2.X = XMIN + (I1 + 0.5) * cellwidth pNewPoint2.Y = YMAX - (I2 - 0.5) * cellheight pNewPoint3.X = XMIN + (I1 - 0.5) * cellwidth pNewPoint3.Y = YMAX - (I2 + 0.5) * cellheight pNewPoint4.X = XMIN + (I1 + 0.5) * cellwidth pNewPoint4.Y = YMAX - (I2 + 0.5) * cellheight 'Get RasterIdentifyObject on these four points Set pIDArray1 = pIdentify.Identify (pNewPoint1) Set pIDArray2 = pIdentify.Identify (pNewPoint2) Set pIDArray3 = pIdentify.Identify (pNewPoint3) 431 Set pIDArray4 = pIdentify.Identify (pNewPoint4) 'Eletotal and Elecount will be used to deal with no surrounding cell or surrounding cell with no data cases Dim Eletotal As Double Eletotal = 0 Dim Elecount As Integer Elecount = 0 'In case one or more surrounding cells not exist, an average of the existing ones is taken If (pIDArray1 Is Nothing) Or (pIDArray2 Is Nothing) Or (pIDArray3 Is Nothing) Or (pIDArray4 Is Nothing) Then If Not pIDArray1 Is Nothing Then Set pRIDObj1 = pIDArray1.Element (0) If pRIDObj1.Name "NoData" Then Elecount = Elecount + 1 Eletotal = Eletotal + CDbl (pRIDObj1.Name) End If End If If Not pIDArray2 Is Nothing Then Set pRIDObj2 = pIDArray2.Element (0) If pRIDObj2.Name "NoData" Then Elecount = Elecount + 1 Eletotal = Eletotal + CDbl (pRIDObj2.Name) End If End If If Not pIDArray3 Is Nothing Then Set pRIDObj3 = pIDArray3.Element (0) If pRIDObj3.Name "NoData" Then Elecount = Elecount + 1 Eletotal = Eletotal + CDbl (pRIDObj3.Name) End If End If If Not pIDArray4 Is Nothing Then Set pRIDObj4 = pIDArray4.Element (0) If pRIDObj4.Name "NoData" Then Elecount = Elecount + 1 Eletotal = Eletotal + CDbl(pRIDObj4.Name) End If End If If Elecount > 0 Then myrow.value (FieldIndex) = Eletotal / Elecount myrow.value (6) = Eletotal / Elecount myrow.Store End If Else Set pRIDObj1 = pIDArray1.Element (0) Set pRIDObj2 = pIDArray2.Element (0) Set pRIDObj3 = pIDArray3.Element (0) Set pRIDObj4 = pIDArray4.Element (0) 'In case one or more surrounding cells have no data, an average is taken 432 If pRIDObj1.Name = "NoData" Or pRIDObj2.Name = "NoData" Or pRIDObj3.Name = "NoData" Or pRIDObj4.Name = "NoData" Then If pRIDObj1.Name "NoData" Then Elecount = Elecount + 1 Eletotal = Eletotal + CDbl (pRIDObj1.Name) End If If pRIDObj2.Name "NoData" Then Elecount = Elecount + 1 Eletotal = Eletotal + CDbl (pRIDObj2.Name) End If If pRIDObj3.Name "NoData" Then Elecount = Elecount + 1 Eletotal = Eletotal + CDbl (pRIDObj3.Name) End If If pRIDObj4.Name "NoData" Then Elecount = Elecount + 1 Eletotal = Eletotal + CDbl (pRIDObj4.Name) End If If Elecount > 0 Then myrow.value (FieldIndex) = Eletotal / Elecount myrow.value (6) = Eletotal / Elecount myrow.Store End If Else 'All 4 surrounding cells exist and have values, elevation for the target point is obtained via bilinear 'interpolation 'E1 to E4 will be used to store the elevations for the surrounding cells Dim E1 As Double Dim E2 As Double Dim E3 As Double Dim E4 As Double 'EE1 and EE2 will be used to store the elevations for the intermediate points Dim EE1 As Double Dim EE2 As Double 'EE will be used to store the elevation for the target point after billinear interpolation Dim EE As Double E1 = CDbl (pRIDObj1.Name) E2 = CDbl (pRIDObj2.Name) E3 = CDbl (pRIDObj3.Name) E4 = CDbl (pRIDObj4.Name) EE1 = E1 + (pNewPoint.X - pNewPoint1.X) * (E2 - E1) / cellwidth EE2 = E3 + (pNewPoint.X - pNewPoint3.X) * (E4 - E3) / cellwidth EE = EE2 + (EE1 - EE2) * (pNewPoint.Y - pNewPoint4.Y) / cellheight myrow.value (FieldIndex) = EE myrow.value (6) = EE myrow.Store End If End If myrow.Store Next j Next i End Sub 433 APPENDIX D PROGRAM CODE FOR CALCULATIONS OF AVERAGE, WEIGHTED AVERAGE, AVERAGE SLOPE CHANGE, AND WEIGHTED AVERAGE SLOPE CHANGE 'Program Name: Slope 'Computes the average slope, weighted average slope, average slope change, and weighted average slope 'change from the 3-D points representing FTSegs in a 3-D space. 'This program was developed by Hubo Cai during th summer of 2003. Option Explicit Sub slope () 'Create a workspace to open the attribute table of the 3-D road centerline data and the target table to store the 'calculations Dim pworkspacefactory As IWorkspaceFactory Set pworkspacefactory = New ShapefileWorkspaceFactory Dim pfeatureworkspace As IFeatureWorkspace Set pfeatureworkspace = pworkspacefactory.OpenFromFile ("c:\dotdata\studyscope", 0) 'Open the 3-D road centerline data (3-D point data) Dim mytable As ITable Set mytable = pfeatureworkspace.OpenTable ("snap_johnston_nc_4ft") MsgBox "how many rows: " & mytable.RowCount (Nothing) 'Open the target table to store the calculations Dim mytargettable As ITable Set mytargettable = pfeatureworkspace.OpenTable ("slope_johnston_nc_average") 'Define all variables Dim sumslopedifference As Double Dim sumslope As Double Dim absoluteslope As Double Dim currentslopedifference As Double Dim slope As Double Dim currentslope As Double Dim current2d As Double Dim currentelevation As Double Dim MERGE As Integer Dim CONDITION As String Dim pointcount As Long Dim segcount As Long Dim sdifference As Double Dim DIST As Double Dim SUM As Double Dim slope1 As Double Dim slope2 As Double Dim SUM_Slope_W As Double Dim SUM_Slope_CH_W As Double 'Loop through the table based on the MERGE value (the ID of the FTSegs) For MERGE = 0 To 100 'Initialize the variables sumslopedifference = 0 sumslope = 0 currentslope = 0 current2d = 0 currentelevation = 0 434 pointcount = 0 DIST = 0 SUM = 0 slope1 = 0 slope2 = 0 SUM_Slope_W = 0 SUM_Slope_CH_W = 0 'Define the tablesort where clause CONDITION = "MERGE = " & MERGE 'Create a tablesort to retrieve all 3-D points belonging to the current FTSeg Dim pTableSort As ITableSort Set pTableSort = New esriCore.TableSort Dim pQueryFilter As IQueryFilter Set pQueryFilter = New QueryFilter pQueryFilter.WhereClause = CONDITION With pTableSort .Fields = "D_TO_S" .Ascending ("D_TO_S") = True Set .QueryFilter = pQueryFilter Set .Table = mytable End With pTableSort.Sort Nothing Dim pCursor As ICursor Set pCursor = pTableSort.Rows Dim pRow As IRow Dim tRow As IRow Set pRow = pCursor.NextRow 'Start from the start point, go through neighboring points and do calculations If Not pRow Is Nothing Then currentelevation = pRow.Value (3) current2d = pRow.Value (5) pointcount = pointcount + 1 Set tRow = mytargettable.CreateRow tRow.Value (1) = pRow.Value (7) Set pRow = pCursor.NextRow Do While Not pRow Is Nothing 'Variable sdifference is the current 2D distance between two neighboring 3D points sdifference = pRow.Value (5) - current2d If sdifference > 0 Then 'DIST is the 3-D distance between two neighboring points DIST = Sqr((pRow.Value (3) - currentelevation) * (pRow.Value (3) - currentelevation) + (pRow.Value (5) - current2d) * (pRow.Value (5) - current2d)) 'SUM is the sum of the 3D distances SUM = SUM + DIST 'Variable slope is the current slope between two neighboring points slope = (pRow.Value (3) - currentelevation) / (pRow.Value (5) - current2d) 'Variable currentslopedifference is the slope change from the current slope to the slope of the previous one currentslopedifference = Abs (slope - currentslope) absoluteslope = Abs (slope) 435 'Variable sumslopedifference holds the sum of the slope changes sumslopedifference = sumslopedifference + currentslopedifference 'Variable sumslope holds the sum of the absolute slopes sumslope = sumslope + absoluteslope 'SUM_Slope_W and SUM_Slope_CH_W hold the sum of the weighted slopes and the sum of the weighted slope changes, respectively SUM_Slope_W = SUM_Slope_W + absoluteslope * DIST SUM_Slope_CH_W = SUM_Slope_CH_W + currentslopedifference * DIST 'MsgBox "2D Difference is: " & (pRow.Value (2) - current2d) & "3D Distance is: " & sum pointcount = pointcount + 1 currentslope = slope Else 'Two points are too close and do nothing MsgBox "Points too close!" End If current2d = pRow.Value (5) currentelevation = pRow.Value (3) Set pRow = pCursor.NextRow Loop 'The totoal number of segments is calculated as the number of 3-D points subtracted by 1 segcount = pointcount - 1 'Calculate all the average values and store them to the target table tRow.Value (2) = sumslopedifference tRow.Value (3) = segcount tRow.Value (4) = sumslopedifference / segcount tRow.Value (5) = sumslope tRow.Value (6) = sumslope / segcount tRow.Value (7) = SUM_Slope_W / SUM tRow.Value (8) = SUM_Slope_CH_W / SUM tRow.Value (9) = SUM tRow.Value (10) = SUM_Slope_W tRow.Value (11) = SUM_Slope_CH_W tRow.Store Else 'No FTSeg with the current MERGE value exists. Do nothing and go to the next MERGE value End If Next End Sub 436 APPENDIX E PROGRAM CODE FOR CHECKING THE ERROR TOLERANCE AND THE MAXIMUM DROP AND MAKING ADJUSTMENTS 'Program Name: Clean_methods 'For all the points along river polylines, check the error tolerance and the maximum drop and make 'corresponding adjustments. 'This program was developed by Hubo Cai during the summer of 2003. Option Explicit Sub Clean_methods() 'Open a workspace and open the point data Dim pworkspacefactory As IWorkspaceFactory Set pworkspacefactory = New ShapefileWorkspaceFactory Dim pfeatureworkspace As IFeatureWorkspace Set pfeatureworkspace = pworkspacefactory.OpenFromFile ("S:\PALRS\users\Hubo\Hydro", 0) Dim mytable As ITable Set mytable = pfeatureworkspace.OpenTable ("PointAlongRiverClean_method1") 'Define variables. Variable totalchanged is used to keep track of the number of points being changed Dim RiverID As Integer Dim elev1 As Double Dim elev2 As Double Dim CONDITION As String Dim totalchanged As Long totalchanged = 0 'Loop through each river For RiverID = 0 To 67 CONDITION = "RiverID = " & RiverID 'elev1 and elev2 are used to store the average elevations of the first 3 points and the last 3 points, respectively elev1 = 0 elev2 = 0 'Define a tablesort Dim pTableSort As ITableSort Set pTableSort = New esriCore.TableSort Dim pQueryFilter As IQueryFilter Set pQueryFilter = New QueryFilter pQueryFilter.whereClause = CONDITION 'Sort 1 to get the first three points With pTableSort .Fields = "D_TO_S" .Ascending ("D_TO_S") = True Set .QueryFilter = pQueryFilter Set .Table = mytable End With pTableSort.Sort Nothing Dim pCursor As ICursor Set pCursor = pTableSort.Rows Dim pRow As IRow Dim trow As IRow Set pRow = pCursor.NextRow 437 Dim count As Integer count = 0 Dim i As Integer 'get the average elevation of the first three points If Not pRow Is Nothing Then Do While (Not pRow Is Nothing) And (count < 3) elev1 = elev1 + pRow.value (1) count = count + 1 Set pRow = pCursor.NextRow Loop elev1 = elev1 / count Else elev1 = 9999 End If 'Sort 2 to get the last three points With pTableSort .Fields = "D_TO_S" .Ascending ("D_TO_S") = False Set .QueryFilter = pQueryFilter Set .Table = mytable End With pTableSort.Sort Nothing Set pCursor = pTableSort.Rows Set pRow = pCursor.NextRow count = 0 'Get the average elevation of the last three points If Not pRow Is Nothing Then Do While (Not pRow Is Nothing) And (count < 3) elev2 = elev2 + pRow.value (1) count = count + 1 Set pRow = pCursor.NextRow Loop elev2 = elev2 / count Else elev2 = 9999 End If Dim down As Boolean Dim upper As Boolean Dim flat As Boolean 'Determine if the water is running from start to the end or from the end to the start If elev1 = 9999 Or elev2 = 999 Then down = False upper = False MsgBox "some points are missing" GoTo eex Else If Abs(elev1 - elev2) <= 10 Then 'difference is less than 1ft down = False upper = False flat = True 438 ElseIf elev1 >= elev2 Then down = True upper = False Else down = False upper = True End If End If 'Start cleaning 'Based on the flowing direction, sort the points based on their distance to the start point so that the order is always from higher elevations to lower elevations If upper Then With pTableSort .Fields = "D_TO_S" .Ascending ("D_TO_S") = False Set .QueryFilter = pQueryFilter Set .Table = mytable End With Else With pTableSort .Fields = "D_TO_S" .Ascending ("D_TO_S") = True Set .QueryFilter = pQueryFilter Set .Table = mytable End With End If pTableSort.Sort Nothing Set pCursor = pTableSort.Rows Set pRow = pCursor.NextRow Dim targetrow As IRow Dim maxdrop As Double Dim tolerance As Double 'Define the maximun drop and the error tolerance maxdrop = 0.1 '10% tolerance = 0.02 '1ft/50ft If Not pRow Is Nothing Then Set targetrow = pRow Set pRow = pCursor.NextRow Do While Not pRow Is Nothing 'check with error tolerance If (pRow.value (1) - targetrow.value (1)) >= 0 Then If Abs((pRow.value (1) - targetrow.value (1)) / ((pRow.value (2) - targetrow.value (2)) * 10)) <= tolerance Then 'Within error tolerance, ok, do nothing Else 'Error tolerance is exceeded. Assign the elevation of the succeeding point to the current point pRow.value (1) = targetrow.value (1) pRow.Store totalchanged = totalchanged + 1 End If Set targetrow = pRow Set pRow = pCursor.NextRow Else 439 'check with maximum drop If Abs((pRow.value (1) - targetrow.value (1)) / ((pRow.value (2) - targetrow.value (2)) * 10)) <= maxdrop Then 'Within the maximum drop range, ok, do nothing Else 'Algorithm 1, adjust the elevation of the current point pRow.value (1) = targetrow.value (1) - maxdrop * Abs(pRow.value (2) - targetrow.value (2)) * 10 pRow.Store totalchanged = totalchanged + 1 'Algorithm 2, adjust the elevation of the succeeding point '(targetrow.value (1) = pRow.value (1) + maxdrop * Abs(pRow.value (2) - targetrow.value (2)) * 10} '{targetrow.Store} '{totalchanged = totalchanged + 1} End If Set targetrow = pRow Set pRow = pCursor.NextRow End If Loop Else End If eex: MsgBox "Some points are missing, pay attention to river with ID: " & RiverID Next End Sub 440 APPENDIX F PROGRAM CODE FOR FLOOD EXTENT PREDICTION 'Program Name: FloodExtentPrediction 'This program is developed to predict flood extent with given LIDAR 20ft DEM and point data along 'waterbodies '(polylines) 'This program was developed by Hubo Cai during the summer of 2003 Option Explicit Sub FloodExtentPrediction () 'Define raster data set and open raster data by calling the OpenRasterDataset Sub Dim praster As IRasterDataset Dim ppraster As IRaster Set ppraster = OpenRasterDataset ("s:\palrs\users\Hubo\Hydro", "elev_wilson").CreateDefaultRaster ppraster.ResampleMethod = RSP_BilinearInterpolation 'Define workspace and open the waterbody data and the resulting polygon data Dim pworkspacefactory As IWorkspaceFactory Set pworkspacefactory = New ShapefileWorkspaceFactory Dim pfeatureworkspace As IFeatureWorkspace Set pfeatureworkspace = pworkspacefactory.OpenFromFile ("S:\PALRS\users\Hubo\Hydro", 0) Dim polygonfeatureclass As IFeatureClass Set polygonfeatureclass = pfeatureworkspace.OpenFeatureClass ("Testing_results_clean_method2") Dim newfeature As IFeature Dim mypointcollection As IPointCollection Dim mypolygon As IPolygon 'a Linefeature represents a waterbody polyline 'linetotal and linecount are used to control the loop Dim linefeatureclass As IFeatureClass Set linefeatureclass = pfeatureworkspace.OpenFeatureClass ("Testing_riv_3") Dim linetotal As Integer linetotal = linefeatureclass.FeatureCount (Nothing) Dim linecount As Integer Dim myfeature As IFeature Dim mypointtable As ITable Set mypointtable = pfeatureworkspace.OpenTable ("PointAlongRiverClean_method2") Dim flood1 As Double Dim flood2 As Double Dim reachingboundary As Long reachingboundary = 0 For linecount = 0 To linetotal - 1 Set myfeature = linefeatureclass.GetFeature (linecount) Dim myline As IPolyline Set myline = myfeature.Shape Dim elevadd As Double 'Assume the given flood is 2 ft 'Our elevation data were multiplied by 10 and therefore, 2 ft equals 20 now elevadd = 20 Dim count As Integer 'Point1, point2, point3, and point4 will be used to represent the vertices 'of the polygon to be constructed. 441 'midpoint represents a point along the waterbody polyline, through which 'the normal line will be built. Dim point1 As IPoint Dim point2 As IPoint Dim point3 As IPoint Dim point4 As IPoint Set point1 = New Point Set point2 = New Point Set point3 = New Point Set point4 = New Point Dim midpoint As IPoint Dim CONDITION As String CONDITION = "RiverID = " & myfeature.value(0) Dim pTableSort As ITableSort Set pTableSort = New esriCore.TableSort Dim pQueryFilter As IQueryFilter Set pQueryFilter = New QueryFilter pQueryFilter.whereClause = CONDITION With pTableSort .Fields = "D_TO_S" .Ascending ("D_TO_S") = True Set .QueryFilter = pQueryFilter Set .Table = mypointtable End With pTableSort.Sort Nothing Dim pCursor As ICursor Set pCursor = pTableSort.Rows Dim pRow As IRow Set pRow = pCursor.NextRow Dim delete1 As Boolean Dim delete2 As Boolean If Not pRow Is Nothing Then 'Built the normal line through the start point of a waterbody polyline feature Dim pnormal1 As ILine Set pnormal1 = New Line myline.QueryNormal esriNoExtension, pRow.value (2), False, 100, pnormal1 Set point1 = pnormal1.ToPoint Dim pnormal2 As ILine Set pnormal2 = New Line Set midpoint = New Point Set midpoint = myline.FromPoint 'The normal line starts from the midpoint to one direction. 'Since we need the normal line in both directions, point2 is 'temporarily used to extent the line to the other direction. point2.X = midpoint.X * 2 - point1.X point2.Y = midpoint.Y * 2 - point1.Y pnormal2.FromPoint = midpoint pnormal2.ToPoint = point2 'Call the mypoint sub to determine point1 and point2 on the normal lines 'The new point1 and point2 having their elevations equaling the target elevation. 442 Call mypoint(delete1, pRow.value (1), point1, ppraster, pnormal1, midpoint, elevadd, flood1) Call mypoint(delete2, pRow.value (1), point2, ppraster, pnormal2, midpoint, elevadd, flood1) If delete1 Then reachingboundary = reachingboundary + 1 Else End If If delete2 Then reachingboundary = reachingboundary + 1 Else End If Dim pnormal3 As ILine Dim pnormal4 As ILine Set pnormal3 = New Line Set pnormal4 = New Line 'Continue to construct normal line and to determine point3 and point4 'at the next point along the water body polyline Set pRow = pCursor.NextRow Do While Not pRow Is Nothing myline.QueryNormal esriNoExtension, pRow.value (2), False, 100, pnormal3 myline.QueryPoint esriNoExtension, pRow.value (2), False, midpoint Set point4 = pnormal3.ToPoint Set pnormal4 = New Line Set point3 = New Point point3.X = midpoint.X * 2 - point4.X point3.Y = midpoint.Y * 2 - point4.Y pnormal4.FromPoint = midpoint pnormal4.ToPoint = point3 'Call the mypoint sub to determine point3 and point4 on the normal lines 'The new point1 and point2 having their elevations equaling the target elevation. Call mypoint(delete1, pRow.value (1), point3, ppraster, pnormal4, midpoint, elevadd, flood2) Call mypoint(delete2, pRow.value (1), point4, ppraster, pnormal3, midpoint, elevadd, flood2) If delete1 Then reachingboundary = reachingboundary + 1 Else End If If delete2 Then reachingboundary = reachingboundary + 1 Else End If 'After 4 vertices are determined, a new polygon is constructed Set mypolygon = New Polygon Set mypointcollection = mypolygon If mypointcollection.pointcount > 0 Then mypointcollection.RemovePoints 0, mypointcollection.pointcount Else End If 443 If delete1 Or delete2 Then If delete1 Then mypointcollection.AddPoints 1, point1 mypointcollection.AddPoints 1, point2 mypointcollection.AddPoints 1, point4 mypolygon.Close Set newfeature = polygonfeatureclass.CreateFeature Set newfeature.Shape = mypolygon newfeature.Store newfeature.value (3) = myfeature.value (0) newfeature.value (4) = (flood1 * 2 + flood2) / 3 newfeature.Store Set point1 = New Point 'Set point2 = New point Set point1 = point4 'Set point2 = point3 flood1 = flood2 Set pRow = pCursor.NextRow Else mypointcollection.AddPoints 1, point1 mypointcollection.AddPoints 1, point2 mypointcollection.AddPoints 1, point3 mypolygon.Close Set newfeature = polygonfeatureclass.CreateFeature Set newfeature.Shape = mypolygon newfeature.Store newfeature.value (3) = myfeature.value(0) newfeature.value (4) = (flood1 * 2 + flood2) / 3 newfeature.Store 'Set point1 = New point Set point2 = New Point 'Set point1 = point4 Set point2 = point3 flood1 = flood2 Set pRow = pCursor.NextRow End If Else mypointcollection.AddPoints 1, point1 mypointcollection.AddPoints 1, point2 mypointcollection.AddPoints 1, point3 mypointcollection.AddPoints 1, point4 mypolygon.Close Set newfeature = polygonfeatureclass.CreateFeature Set newfeature.Shape = mypolygon newfeature.Store newfeature.value (3) = myfeature.value (0) newfeature.value (4) = (flood1 + flood2) / 2 newfeature.Store 'After the new polygon is constructed and stored, point3 and point4 'are assigned to point2 and point1, respectively. 'The procedure continues to identify new point3 and point4 at the next 'point along the water body polyline and consequently, new polygon will 'be constructed. 444 'Doing the inner loop will construct polygons surrounding the water body. Set point1 = New Point Set point2 = New Point Set point1 = point4 Set point2 = point3 flood1 = flood2 Set pRow = pCursor.NextRow End If Loop Else End If Next MsgBox "total reaching boundary: " & reachingboundary End Sub Public Function OpenRasterDataset (sPath As String, sFileName As String) As IRasterDataset ' sPath: directory where dataset resides ' sFileName: name of the raster dataset On Error GoTo ErrorHandler ' Create RasterWorkSpaceFactory Dim pWSF As IWorkspaceFactory Set pWSF = New RasterWorkspaceFactory ' Get RasterWorkspace Dim pRasWS As IRasterWorkspace If pWSF.IsWorkspace (sPath) Then Set pRasWS = pWSF.OpenFromFile (sPath, 0) Set OpenRasterDataset = pRasWS.OpenRasterDataset (sFileName) End If ' Release memeory Set pRasWS = Nothing Set pWSF = Nothing Exit Function ErrorHandler: Set OpenRasterDataset = Nothing End Function 'This sub aims at determining two points on the normal line, which have their elevations 'equalling to the target elevation. Public Sub mypoint(delete As Boolean, base_elev As Double, outpoint As IPoint, pInRaster As IRaster, pInLine As ILine, ppInPoint As IPoint, elevadd As Double, flood As Double) ' pInRaster: input raster ' pInFeatureClass: input line feature class ' sFieldName: name of the field that stores the values delete = False ' Create a raster layer and QI for IIdentify interface Dim pRLayer As IRasterLayer Set pRLayer = New RasterLayer pInRaster.ResampleMethod = RSP_BilinearInterpolation pRLayer.CreateFromRaster pInRaster Dim pIdentify As IIdentify Set pIdentify = pRLayer Dim pProperty As IRasterProps Set pProperty = pInRaster Dim cellheight As Double cellheight = pProperty.Height 445 cellheight = (pProperty.Extent.YMAX - pProperty.Extent.YMin) / cellheight 'Determine cell width Dim cellwidth As Double cellwidth = pProperty.Width cellwidth = (pProperty.Extent.XMax - pProperty.Extent.XMIN) / cellwidth 'Determine the upper-left corner of the raster data Dim XMIN As Double XMIN = pProperty.Extent.XMIN Dim YMAX As Double YMAX = pProperty.Extent.YMAX 'I1 and I2 will be used to locate the surrounding cells for the target point 'and determine the centerpoints for these surrounding cells Dim I1 As Integer Dim I2 As Integer Dim D1 As Double Dim D2 As Double 'pIDArray and pRIDObj are used for the target point while the others are used for the surrounding cell centers. Dim pIDArray As IArray Dim pRIDObj As IRasterIdentifyObj Dim pIDArray1 As IArray Dim pRIDObj1 As IRasterIdentifyObj Dim pIDArray2 As IArray Dim pRIDObj2 As IRasterIdentifyObj Dim pNewPoint1 As IPoint Set pNewPoint1 = New Point Dim pNewPoint2 As IPoint Set pNewPoint2 = New Point 'Get RasterIdentifyObject on the point on the river line Dim pInPoint As IPoint Set pInPoint = New Point Set pInPoint = ppInPoint 'MsgBox "Do we have pInPoint? " & pInPoint.X & ", " & pInPoint.Y Set pIDArray = pIdentify.Identify(pInPoint) If pIDArray Is Nothing Then MsgBox "No elevation is found, exit sub!" Exit Sub Else 'MsgBox "GOOD POINT!" End If 'Determine the elevation of the midpoint and consequently 'the target elevation 'Set pRIDObj = pIDArray.Element (0) 'Dim base_elev As Double 'base_elev = CDbl (pRIDObj.Name) Dim target_elev As Double target_elev = base_elev + elevadd Dim currentelev As Double Dim elev1 As Double Dim elev2 As Double 446 'Identify two points on line with one elevation lower than the target elevation 'and the other with elevation higher than the target elevation Dim ii As Integer ii = 0 'The Do loop searches for the target point along the normal line, starting from 'the mid point. Each search step has the distance increase equalling to cell 'width. The target point is located via linear interpolation. Do ii = ii + 1 pInLine.QueryPoint esriExtendTangentAtTo, ii * cellwidth, False, pNewPoint2 Set pIDArray2 = pIdentify.Identify (pNewPoint2) If pIDArray2 Is Nothing Then 'MsgBox "ii is: " & ii Exit Do Else Set pRIDObj2 = pIDArray2.Element (0) If pRIDObj2.Name = "NoData" Then 'MsgBox "new ii is: " & ii Exit Do Else currentelev = CDbl (pRIDObj2.Name) elev2 = currentelev End If End If Loop While (currentelev "NoData") pInLine.QueryPoint esriExtendTangentAtTo, (ii - 1) * cellwidth, False, pNewPoint1 Set pIDArray1 = pIdentify.Identify (pNewPoint1) Set pRIDObj1 = pIDArray1.Element (0) elev1 = CDbl (pRIDObj1.Name) Set outpoint = New Point If pIDArray2 Is Nothing Or pRIDObj2.Name = "NoData" Then 'The boundary is reached delete = True outpoint.X = pNewPoint1.X outpoint.Y = pNewPoint1.Y Else delete = False If Abs(elev2 - elev1) < 0.01 Then 'Two points having elevation very close to each other 'the target point is located in the middle of the 'two points outpoint.X = (pNewPoint1.X + pNewPoint2.X) / 2 outpoint.Y = (pNewPoint1.Y + pNewPoint2.Y) / 2 Else 'Linear interpolation is applied outpoint.X = pNewPoint1.X + (pNewPoint2.X - pNewPoint1.X) * (target_elev - elev1) / (elev2 - elev1) outpoint.Y = pNewPoint1.Y + (pNewPoint2.Y - pNewPoint1.Y) * (target_elev - elev1) / (elev2 - elev1) End If End If 'This flood is the flooded water level of the new polygon being constructed. flood = target_elev 'MsgBox "Base elevatio is: " & base_elev & "; flood elevation is: " & target_elev End Sub 447 APPENDIX G PROGRAM CODE FOR FLOODE ROAD SEGMENT IDENTIFICATION 'Program Name: FloodedRoadSegmentIdentification 'This program is developed to identify flooded road segments given the resulting polygons representing flood 'extent 'and the 3-D road data. 'This program was developed by Hubo Cai during the summer of 2003 Option Explicit 'The main sub determines the road segments in the flood extents. For each of the road segments, the public sub 'myflood is called to determine the flooded road segments within the road segments in the flood extent. Sub FloodedRoadSegmentIdentification () 'Open workspace and all the files Dim pworkspacefactory As IWorkspaceFactory Set pworkspacefactory = New ShapefileWorkspaceFactory Dim pfeatureworkspace As IFeatureWorkspace Set pfeatureworkspace = pworkspacefactory.OpenFromFile ("s:\palrs\users\hubo\Hydro", 0) 'Open 3-D road data file Dim linefeatureclass As IFeatureClass Set linefeatureclass = pfeatureworkspace.OpenFeatureClass ("Inter_Wilson") 'Open the flood extent file Dim polygonfeatureclass As IFeatureClass Set polygonfeatureclass = pfeatureworkspace.OpenFeatureClass ("Testing_results_500_4") 'Open the new road segment file to store identified flooded road segments Dim newlinefeatureclass As IFeatureClass Set newlinefeatureclass = pfeatureworkspace.OpenFeatureClass ("road1") Dim mylinefeature As IFeature Dim myline As IPolyline 'Variable linecount is used to identify the number of polylines in the 3-D road data Dim linecount As Long linecount = linefeatureclass.FeatureCount (Nothing) Dim flood As Double Dim lineID As Long 'The outer loop works with each feature from the 3-D road data For lineID = 0 To linecount - 1 Set mylinefeature = linefeatureclass.GetFeature (lineID) Set myline = New Polyline Set myline = mylinefeature.Shape Dim myoperator As ITopologicalOperator Set myoperator = myline Dim mypolygonfeature As IFeature 'The variable polygoncount is used to count the number of polygons in the flood extent file Dim polygoncount As Long polygoncount = polygonfeatureclass.FeatureCount(Nothing) Dim polygonID As Long 'The inner loop works with every polygon from the flood extent file For polygonID = 0 To polygoncount - 1 Set mypolygonfeature = polygonfeatureclass.GetFeature (polygonID) 'Obtain the flood level stored with each polygon flood = mypolygonfeature.value (4) / 10 Dim mypolygon As IPolygon 448 Set mypolygon = New Polygon Set mypolygon = mypolygonfeature.Shape Dim newline As IPolyline Set newline = New Polyline Dim newfeature As IFeature 'Object newmultipoints will be used to store points for construcing flooded road segments Dim newmultipoints As IMultipoint Set newmultipoints = New Multipoint Dim newpointcollection As IPointCollection Dim mygeometry As IGeometry Set mygeometry = New Multipoint 'Each road feature is intersected with each polygon Set mygeometry = myoperator.intersect (mypolygon, 1) Set newpointcollection = mygeometry If newpointcollection.pointcount = 0 Then Else 'Work with the intersecting points to figure out the points 'needed to construct flooded road segments Dim ii As Integer ii = newpointcollection.pointcount Mod 2 If ii = 0 Then Else 'identify start and end nodes 'Judge if start or end nodes are in the polygon 'Determine if the start or end node should be used to construct flooded road segments Dim startnode As IPoint Dim endnode As IPoint Set startnode = New Point Set endnode = New Point Set startnode = myline.FromPoint Set endnode = myline.ToPoint Dim test As IRelationalOperator Set test = mypolygon Dim startin As Boolean Dim endin As Boolean startin = test.Contains (startnode) endin = test.Contains (endnode) If startin Then newpointcollection.AddPoint startnode ElseIf endin Then newpointcollection.AddPoint endnode Else MsgBox "something wrong!" MsgBox "line ID is:" & mylinefeature.value (0) MsgBox "polygon ID is:" & mypolygonfeature.value (0) Exit Sub End If End If Dim pointcount As Integer pointcount = newpointcollection.pointcount 'The two arrays defined hereby will be used to store those points 'in the order of their distances to the start point of the road segment 'from which they are identified. 'By doing this, all intersected road segments are identified. 449 Dim orderedpointcollection (100) As IPoint Dim dist (100) As Double Dim i As Integer Dim out As IPoint Dim dis1 As Double Dim dis2 As Double Dim bright As Boolean For i = 0 To pointcount - 1 Set out = New Point myline.QueryPointAndDistance esriNoExtension, newpointcollection.Point (i), False, out, dis1, dis2, bright dist (i) = dis1 Next Dim outerloop As Integer Dim interloop As Integer Dim min As Double For outerloop = 0 To pointcount - 2 min = dist (outerloop) For interloop = outerloop + 1 To pointcount - 1 If dist (interloop) < min Then min = dist (interloop) dist (interloop) = dist (outerloop) dist (outerloop) = min Else End If Next Next Dim newlinecount As Integer newlinecount = pointcount / 2 Dim linecreation As Integer For linecreation = 0 To newlinecount - 1 Set newline = New Polyline 'Call the myflood function to determine the flooded road segments 'for all road segments that are in the flood extent Call myflood(mylinefeature.value (2), dist(linecreation * 2), dist(linecreation * 2 + 1), flood, myline, newlinefeatureclass) Next End If Next Next End Sub 'This public works with the 3-D road data to determine flooded road segments within the road segments in the 'flood extent Public Sub myflood(linemerge As Long, dist1 As Double, dist2 As Double, flood As Double, myline As IPolyline, newlinefeatureclass As IFeatureClass) 'Open workspace and needed files Dim fpworkspacefactory As IWorkspaceFactory Set fpworkspacefactory = New ShapefileWorkspaceFactory Dim fpfeatureworkspace As IFeatureWorkspace 450 Set fpfeatureworkspace = fpworkspacefactory.OpenFromFile ("S:\PALRS\users\Hubo\Nroad_english\WILSON", 0) 'This table contains 3-D road data Dim fmytable As ITable Set fmytable = fpfeatureworkspace.OpenTable("snap_wilson") Dim newfeature As IFeature 'Variable CONDITION is that condition that will be used to retrieve all 3-D points enclosed 'by the start and end points of a road segment in the flood extent Dim CONDITION As String CONDITION = "MERGE = " & linemerge & " AND D_TO_S >= " & dist1 & " AND D_TO_S <= " & dist2 'Variables CONDITION1 and CONDITION2 will store the necessary conditions for retrieving 3-D 'points from the table, which are used to determine the elevations for the start and end points 'of the road segment in the flood extent, via linear interpolation from two neighboring points. Dim CONDITION1 As String Dim CONDITION2 As String CONDITION1 = "MERGE = " & linemerge & " AND D_TO_S <= " & dist1 CONDITION2 = "MERGE = " & linemerge & " AND D_TO_S >= " & dist2 Dim pTableSort As ITableSort Set pTableSort = New esriCore.TableSort Dim pQueryFilter As IQueryFilter Set pQueryFilter = New QueryFilter pQueryFilter.whereClause = CONDITION1 Dim newline As IPolyline Set newline = New Polyline Dim e1 As Double Dim e2 As Double Dim ss1 As Double Dim ss2 As Double 'Obtain the elevation and distance the 3-D point right before dist1 (the start 'point of the road segment in the flood extent With pTableSort .Fields = "D_TO_S" .Ascending ("D_TO_S") = True Set .QueryFilter = pQueryFilter Set .Table = fmytable End With pTableSort.Sort Nothing Dim pCursor As ICursor Set pCursor = pTableSort.Rows Dim pRow As IRow Set pRow = pCursor.NextRow Do While Not pRow Is Nothing e1 = pRow.value (3) ss1 = pRow.value (5) Set pRow = pCursor.NextRow Loop 'Obtain the elevation and distance to the start point for the point right after dist2 pQueryFilter.whereClause = CONDITION2 With pTableSort 451 .Fields = "D_TO_S" .Ascending ("D_TO_S") = True Set .QueryFilter = pQueryFilter Set .Table = fmytable End With pTableSort.Sort Nothing Set pCursor = Nothing Set pRow = Nothing Set pCursor = pTableSort.Rows Set pRow = pCursor.NextRow If (Not pRow Is Nothing) Then e2 = pRow.value (3) ss2 = pRow.value (5) Else MsgBox "Something wrong!" Exit Sub End If 'Finding 3-D points between dist1 and dist2 pQueryFilter.whereClause = CONDITION With pTableSort .Fields = "D_TO_S" .Ascending ("D_TO_S") = True Set .QueryFilter = pQueryFilter Set .Table = fmytable End With pTableSort.Sort Nothing Set pCursor = Nothing Set pRow = Nothing Set pCursor = pTableSort.Rows Set pRow = pCursor.NextRow Dim EE1 As Double Dim nopoint As Boolean If (Not pRow Is Nothing) Then 'Determining the elevation for the start point via linear interpolation Dim difference1 As Double difference1 = pRow.value (5) - ss1 If Abs (difference1) < 0.001 Then EE1 = pRow.value (3) Else EE1 = e1 + (pRow.value (3) - e1) * (dist1 - ss1) / (pRow.value (5) - ss1) End If Else 'In case there is no point between dist1 and dist2, the whole road segment between 'dist1 and dist2 is assumed to be flooded nopoint = True Set newline = Nothing If dist1 <= dist2 Then myline.GetSubcurve dist1, dist2, False, newline Set newfeature = newlinefeatureclass.CreateFeature Set newfeature.Shape = newline newfeature.value (2) = linemerge newfeature.value (3) = dist1 452 newfeature.value (4) = dist2 newfeature.Store 'MsgBox "create line at position 1" Exit Sub Else End If End If With pTableSort .Fields = "D_TO_S" .Ascending("D_TO_S") = False Set .QueryFilter = pQueryFilter Set .Table = fmytable End With pTableSort.Sort Nothing Set pCursor = Nothing Set pRow = Nothing Set pCursor = pTableSort.Rows Set pRow = pCursor.NextRow Dim EE2 As Double 'Determining the elevation of the end point If (Not pRow Is Nothing) Then Dim difference2 As Double difference2 = ss2 - pRow.value (5) If Abs (difference1) < 0.001 Then EE2 = pRow.value (3) Else EE2 = e2 - (e2 - pRow.value (3)) * (dist2 - ss2) / (ss2 - pRow.value (5)) End If Else 'There is no point between dist1 and dist2 and shoule already be dealt with Exit Sub End If 'Start to identify flooded road segments for the road segment in the flood extent 'substart and subend are used to label the start and end points for the flooded segments Dim substart As Double Dim subend As Double Dim previousrow As IRow Dim currentrow As IRow Dim afterrow As IRow If EE1 < flood Then 'ss1 is the start point for the first flooded segment substart = ss1 Else 'If ss1 is not the start point for the first flooded segment, 'the first point having elevation lower than the flood water level is identified 'the first point having elevation equaling to the flood water level is identified vai 'linear interpolation With pTableSort .Fields = "D_TO_S" .Ascending ("D_TO_S") = True Set .QueryFilter = pQueryFilter 453 Set .Table = fmytable End With pTableSort.Sort Nothing Set pCursor = Nothing Set pRow = Nothing Set pCursor = pTableSort.Rows Set pRow = pCursor.NextRow If Not pRow Is Nothing Then Set currentrow = pRow Else Exit Sub End If Do While (Not pRow Is Nothing) If pRow.value (3) > flood Then Set currentrow = pRow Set pRow = pCursor.NextRow Else Exit Do End If Loop If (Not pRow Is Nothing) Then 'afterrow indicates the first point having its elevation lower than the flood water level Set afterrow = pRow If (pRow Is currentrow) Then 'This indicates the start point is the point right before afterrow substart = ss1 + (afterrow.value (5) - ss1) * (EE1 - flood) / (EE1 - afterrow.value (3)) Else 'This indicates the current row represents the point right before afterrow Set previousrow = currentrow If Abs(previousrow.value (3) - afterrow.value (3)) < 0.0001 Then substart = previousrow.value (5) Else substart = previousrow.value(5) + (afterrow.value(5) - previousrow.value(5)) * (previousrow.value(3) - flood) / (previousrow.value(3) - afterrow.value(3)) End If End If Else 'This indicates the failure to find a point having its elevation lower than the flood water level 'and therefore, the sub is ended. There is no flooded road segment Exit Sub End If End If 'After identifying the sub start point, it continues to search for the sub end point 'tempRow is a variable used to store the current position of pRow Dim tempRow As IRow Set tempRow = pRow Set pRow = pCursor.NextRow If (Not pRow Is Nothing) Then Set currentrow = pRow Else 'Reaches the end 'use the ss2 point to get the flooded road segment 454 If EE2 <= flood Then subend = ss2 Else subend = tempRow.value(5) + (flood - tempRow.value(3)) * (ss2 - tempRow.value(5)) / (EE2 - tempRow.value(3)) End If 'Creating a linear feature for the identified flooded road segment If substart <= subend Then myline.GetSubcurve substart, subend, False, newline Set newfeature = newlinefeatureclass.CreateFeature Set newfeature.Shape = newline newfeature.value (2) = linemerge newfeature.value (3) = substart newfeature.value (4) = subend newfeature.Store Set newline = Nothing Else Exit Sub End If Exit Sub End If 'It continues to find the first point with elevation higher than the flood water level Do While (Not pRow Is Nothing) If pRow.value (3) < flood Then Set currentrow = pRow Set pRow = pCursor.NextRow 'MsgBox "currentrow is not nothing" & currentrow.Value (5) Else Exit Do End If Loop 'Determining the sub end point, which has its elevation equaling to the flood water level If (Not pRow Is Nothing) Then Set afterrow = pRow If (currentrow Is pRow) Then Set previousrow = tempRow Else Set previousrow = currentrow End If 'Set afterrow = pRow If Abs(previousrow.value(3) - afterrow.value(3)) < 0.0001 Then subend = afterrow.value(5) Else subend = previousrow.value(5) + (afterrow.value(5) - previousrow.value(5)) * (previousrow.value(3) - flood) / (previousrow.value(3) - afterrow.value(3)) End If 'Creat a new line feature for the identified flooded road segment If substart <= subend Then myline.GetSubcurve substart, subend, False, newline Set newfeature = newlinefeatureclass.CreateFeature Set newfeature.Shape = newline newfeature.value (2) = linemerge 455 newfeature.value (3) = substart newfeature.value (4) = subend newfeature.Store 'MsgBox "create line at position 2" Set newline = Nothing Else End If Else 'It indicates that it reaches the end, but doesn't find a point with elevation 'higher than the flood water level 'This means that ss2 and currentrow should be used to determine the flooded road segment If EE2 <= flood Then subend = ss2 Else subend = currentrow.value(5) + (flood - currentrow.value(3)) * (ss2 - currentrow.value(5)) / (EE2 - currentrow.value(3)) End If If substart <= subend Then myline.GetSubcurve substart, subend, False, newline Set newfeature = newlinefeatureclass.CreateFeature Set newfeature.Shape = newline newfeature.value (2) = linemerge newfeature.value (3) = substart newfeature.value (4) = subend newfeature.Store Set newline = Nothing Else End If Exit Sub End If 'Continue the loops to find the next flooded road segment Do While (Not pRow Is Nothing) Set currentrow = pRow Set pRow = pCursor.NextRow Do While (Not pRow Is Nothing) If pRow.value(3) > flood Then Set currentrow = pRow Set pRow = pCursor.NextRow Else Exit Do End If Loop If (Not pRow Is Nothing) Then Set previousrow = currentrow Set afterrow = pRow If Abs(previousrow.value(3) - afterrow.value(3)) < 0.0001 Then substart = previousrow.value(5) Else substart = previousrow.value(5) + (afterrow.value(5) - previousrow.value(5)) * (previousrow.value(3) - flood) / (previousrow.value(3) - afterrow.value(3)) End If Set currentrow = pRow 456 Else If EE2 < flood Then subend = ss2 substart = currentrow.value(5) + (flood - currentrow.value(3)) * (ss2 - currentrow.value(5)) / (EE2 - currentrow.value(3)) If substart <= subend Then myline.GetSubcurve substart, subend, False, newline Set newfeature = newlinefeatureclass.CreateFeature Set newfeature.Shape = newline newfeature.value (2) = linemerge newfeature.value (3) = substart newfeature.value (4) = subend newfeature.Store Set newline = Nothing Else End If Else End If Exit Sub End If Set pRow = pCursor.NextRow If pRow Is Nothing Then If ss2 <= flood Then subend = ss2 Else subend = currentrow.value(5) + (ss2 - currentrow.value(5)) * (flood - currentrow.value(3)) / (EE2 - currentrow.value(3)) End If Set newline = Nothing If substart <= subend Then myline.GetSubcurve substart, subend, False, newline Set newfeature = newlinefeatureclass.CreateFeature Set newfeature.Shape = newline newfeature.value (2) = linemerge newfeature.value (3) = substart newfeature.value (4) = subend newfeature.Store Set newline = Nothing Exit Sub Else End If Else Do While (Not pRow Is Nothing) If pRow.value (3) < flood Then Set currentrow = pRow Set pRow = pCursor.NextRow Else Exit Do End If Loop End If If Not pRow Is Nothing Then Set previousrow = currentrow 457 Set afterrow = pRow If Abs(previousrow.value (3) - afterrow.value (3)) < 0.0001 Then subend = afterrow.value(5) Else subend = previousrow.value(5) + (afterrow.value(5) - previousrow.value(5)) * (previousrow.value(3) - flood) / (previousrow.value(3) - afterrow.value(3)) End If Else If EE2 > flood Then subend = currentrow.value(5) + (ss2 - currentrow.value(5)) * (flood - currentrow.value(3)) / (EE2 - currentrow.value(3)) Else Exit Sub End If End If Set newline = Nothing If substart <= subend Then myline.GetSubcurve substart, subend, False, newline Set newfeature = newlinefeatureclass.CreateFeature Set newfeature.Shape = newline newfeature.value (2) = linemerge newfeature.value (3) = substart newfeature.value (4) = subend newfeature.Store Set newline = Nothing Else End If Loop End Sub