The VBA code provided at the bottom of this document is an updated
version (from ArcGIS 9.0 to ArcGIS 9.2) of the polygon smoothing
algorithm described below. A bug that occurred when multiple wells had
the same location was also fixed.

# 

# SMOOTH OIL & GAS FIELD OUTLINE POLYGONS

# MADE FROM BUFFERED WELLS

Why smooth buffered field outlines? See the issues in the figure below:

![](media/image2.wmf)

The smoothing application provided as VBA code below does the following:

-   Adds area to the concave portions; doesn't add area to convex
    portions to maintain buffer spacing

-   Fills in non-field "islands" smaller than buffer size

-   Joins separate polygon rings with a "bridge" if sufficiently close

-   Minimizes increase in total field area

Methodology: creates trapezoids between neighboring wells within an
oil/gas field. Wells are considered neighbors if the distance between
them is less than or equal to 2.5x the average of their radii (see
figure below of "rules").

Creates a new \"smoothed\" polygon feature class (FC) from the second
layer in the map. The new FC is in the same workspace as the original
FC, but with \"\_smu\" suffix added to name.

The first layer must contain the well points originally used to create
the unsmoothed field outlines (which are in the second layer).

The original polygon outline is used to optimize the \"unioning\" step.
If a trapezoid falls entirely within the original field outline then it
is omitted since it will not change the outline.

![](media/image2.wmf)

Results of smoothing application as applied to the unsmoothed above
example:

![](media/image2.wmf)

Copy the code into ThisDocument or a module.

Set up a UIControl button to run the code.

Inputs:

1)  Layer (0) is the well points used to make the field outline

2)  Layer (1) is unsmoothed buffered field outlines

Output: Smoothed outlines with suffix "\_smu" in the same workspace as
unsmoothed outlines.

Code by Kirk Kuykendall, AmberGIS; text by Sam Limerick, Z, Inc.

==============================================================

Option Explicit

Private m_TC As ITrackCancel

Private Declare Function GetTickCount Lib \"kernel32\" () As Long

Sub Smooth()

\'

\' Creates a new \"smoothed\" polygon featureclass from the second

\' layer in the TOC(which is an \"unsmoothed\" polygon)in the map.

\' The new featureclass is created in the same workspace as the

\' original featureclass, but with a \"\_smu\" suffix. You can use

\' shapefiles or GDB feature classes.

\'

\' The first layer in the TOC must contain the points originally used

\' to create the unsmoothed polyon in layer 2(layer 2 created by
buffering

\' the points with a radius and unioning the resultant polygons).

\'

\' The approach involves creating trapezoids between neighboring

\' wells. Wells are considered neighbors if the distance between

\' them is less than or equal to 2.5x the average of their radii.

\'

\' The original polygon outline is used to optimize the \"unioning\"

\' step. If a trapezoid falls entirely within the original field

\' outline then it is omitted since it will not change the outline

\' any.

\'

\' Watch the lower left corner task bar for some progress reports

\' that the processing is indeed happening.

\'

On Error GoTo EH

Dim pMxDoc As IMxDocument

Set pMxDoc = ThisDocument

Dim pMC As IMouseCursor

Set pMC = New MouseCursor

pMC.SetCursor 2

Dim lTicks As Long

lTicks = GetTickCount

\' get the wells point layer

Dim pWellFLayer As IFeatureLayer

Set pWellFLayer = pMxDoc.FocusMap.Layer(0)

If pWellFLayer.FeatureClass.ShapeType \<\> esriGeometryPoint Then

MsgBox \"first layer must be a point layer\"

Exit Sub

End If

\' make sure the required columns are present on the wells featureclass

Dim vFld As Variant

For Each vFld In Array(\"buffer\", \"field\")

If pWellFLayer.FeatureClass.FindField(vFld) = -1 Then

MsgBox vFld & \" column not found on \" & pWellFLayer.Name

Exit Sub

End If

Next vFld

\' get the layer with the unsmoothed field outlines

Dim pBufferLayer As IFeatureLayer

Set pBufferLayer = pMxDoc.FocusMap.Layer(1)

If pBufferLayer.FeatureClass.ShapeType \<\> esriGeometryPolygon Then

MsgBox \"second layer must be a polygon layer\"

Exit Sub

End If

\' make sure the required columns are present on the field outline

\' featureclass.

For Each vFld In Array(\"field\")

If pBufferLayer.FeatureClass.FindField(vFld) = -1 Then

MsgBox vFld & \" column not found on \" & pBufferLayer.Name

Exit Sub

End If

Next vFld

\' if there are selected field outlines, then just process those,

\' otherwise process all features in the unsmoothed layer.

Dim pInFCur As IFeatureCursor

Dim pFSel As IFeatureSelection

Set pFSel = pBufferLayer

If pFSel.SelectionSet.Count = 0 Then

Set pInFCur = pBufferLayer.FeatureClass.Search(Nothing, False)

Else

pFSel.SelectionSet.Search Nothing, False, pInFCur

End If

\' create a new featureclass that has same columns as the unsmoothed

\' featureclass, overwriting it if it already exists.

\' (also adds a \"PctChange\" field)

Dim pOutFC As IFeatureClass

Set pOutFC = CreateFC(pBufferLayer.FeatureClass, \_

pBufferLayer.FeatureClass.AliasName & \"\_smu\")

If pOutFC Is Nothing Then

MsgBox \"no output featureclass created\"

Exit Sub

End If

If pOutFC.FindField(\"PctChange\") = -1 Then

MsgBox \"new featureclass is missing PctChange field\"

Exit Sub

End If

Set m_TC = New CancelTracker

m_TC.CancelOnKeyPress = True

Dim pOutFCur As IFeatureCursor

Set pOutFCur = pOutFC.Insert(True)

Dim pBufferFeat As IFeature

Set pBufferFeat = pInFCur.NextFeature

Do Until pBufferFeat Is Nothing

\' put the name of the field being processed in the caption, that

\' way if there is a crash, we can select that field and

\' step through the program and debug it more easily.

Application.Caption = pBufferFeat.Value(pInFCur.FindField(\"field\"))

\' make a smoothed field outline using the unsmoothed outline

\' plus the wells belonging to the oilfield.

Dim pFieldOutline As IPolygon

Set pFieldOutline = MakeFieldOutline(pWellFLayer.FeatureClass, \_

pBufferFeat)

If Not m_TC.Continue Then Exit Do

If Not pFieldOutline Is Nothing Then

\' compute percent change

Dim dPctChange As Double

dPctChange = GetArea(pFieldOutline) / GetArea(pBufferFeat.Shape)

Dim pFeatBuff As IFeatureBuffer

Set pFeatBuff = pOutFC.CreateFeatureBuffer

Set pFeatBuff.Shape = pFieldOutline

\' copy all column values

CopyColumns pFeatBuff, pBufferFeat

\' moved to avoid problems if pctchange is present in original featclass

pFeatBuff.Value(pFeatBuff.Fields.FindField(\"PctChange\")) = dPctChange

pOutFCur.InsertFeature pFeatBuff

Else

Debug.Print \"nothing written for field \" & \_

pBufferFeat.Value(pBufferFeat.Fields.FindField(\"Field\"))

End If

Set pBufferFeat = pInFCur.NextFeature

Loop

pOutFCur.Flush

If m_TC.Continue Then

\' wake up the user

Beep

Application.Caption = \"\"

lTicks = GetTickCount - lTicks

MsgBox \"Done seconds = \" & lTicks / 1000

Else

MsgBox \"canceled\"

End If

Set m_TC = Nothing

Exit Sub

EH:

MsgBox \"error in Smooth: \" & Err.Description & \", \" & Err.Number

Set m_TC = Nothing

End Sub

Function MakeFieldOutline(pWellFC As IFeatureClass, \_

pBufferFeat As IFeature) As IPolygon

On Error GoTo EH

Dim sField As String

sField = pBufferFeat.Value(pBufferFeat.Fields.FindField(\"Field\"))

Dim sWhere As String

sWhere = \"field = \'\" & sField & \"\'\"

Dim cSubFieldIDs As Collection

Set cSubFieldIDs = GetSubFields(pWellFC, sWhere)

Dim l As Long

Dim pOutPolygon As IPolygon

For l = 1 To cSubFieldIDs.Count

Dim pPolygon As IPolygon

Set pPolygon = MakeSubFieldOutline(pWellFC, pBufferFeat, sField,
cSubFieldIDs.Item(l))

If pOutPolygon Is Nothing Then

Set pOutPolygon = pPolygon

Else

Dim pTopoOp As ITopologicalOperator

Set pTopoOp = pOutPolygon

Set pOutPolygon = pTopoOp.Union(pPolygon)

End If

If Not m_TC.Continue Then Exit Function

Next l

Set MakeFieldOutline = pOutPolygon

Exit Function

EH:

MsgBox \"error in MakeFieldOutlines: \" & Err.Description & \", \" &
Err.Number

End Function

Function GetSubFields(pWellFC As IFeatureClass, sWhere As String) As
Collection

\'

\' return a collection of unique subField ID values.

\'

Dim pOutColl As Collection

Set pOutColl = New Collection

Dim lFld As Long

lFld = pWellFC.FindField(\"SubField\")

If lFld = -1 Then

\' no subfield column so just add one id (0)

pOutColl.Add -1

Set GetSubFields = pOutColl

Exit Function

End If

Dim pQF As IQueryFilter

Set pQF = New QueryFilter

pQF.WhereClause = sWhere

Dim pFCur As IFeatureCursor

Set pFCur = pWellFC.Search(pQF, False)

Dim pDict As Scripting.Dictionary

Set pDict = New Dictionary

Dim pFeat As IFeature

Set pFeat = pFCur.NextFeature

Do Until pFeat Is Nothing

Dim lSubFieldID As Long

lSubFieldID = pFeat.Value(lFld)

If pDict.Exists(lSubFieldID) Then

pDict.Item(lSubFieldID) = pDict.Item(lSubFieldID) + 1

Else

pDict.Add lSubFieldID, 1

End If

Set pFeat = pFCur.NextFeature

Loop

Dim vID As Variant

For Each vID In pDict.Keys

pOutColl.Add vID

Next

Set GetSubFields = pOutColl

End Function

Function MakeSubFieldOutline(pWellFC As IFeatureClass, \_

pBufferFeat As IFeature, \_

ByVal sOilFieldName As String, \_

ByVal lSubFieldID As Long) As IPolygon

\'

\' union more trapezoids with pBufferFeat.Shape to and

\' return the unioned polygon. Such that the trapezoids are formed

\' by pairs of neighboring wells.

\' find max buffer radius

On Error GoTo EH

Dim sWhere As String

If lSubFieldID \> -1 Then

sWhere = \"field = \'\" & sOilFieldName & \"\' and Subfield = \" &
lSubFieldID

Else

sWhere = \"field = \'\" & sOilFieldName & \"\'\"

End If

Debug.Print \"finding all wells \" & sWhere

\' make a hashtable of radii keyed by OID, and maxradius

Dim dMaxRadius As Double

Dim pRadiusDict As Scripting.Dictionary

GetRadii pWellFC, sWhere, pRadiusDict, dMaxRadius

If dMaxRadius \< 0 Then Exit Function

If Not m_TC.Continue Then Exit Function

Dim lFld As Long

lFld = pWellFC.FindField(\"buffer\")

Dim pPntDict As Scripting.Dictionary

Set pPntDict = LoadDictionary(pWellFC, sWhere)

\' Its called a matrix since it can be thought of as a symmetric binary

\' sparse matrix.

\' populate pOIDDict with potential neighbors, later we will cull some

\' of the neighbors.

Dim pOIDDict As Scripting.Dictionary

Set pOIDDict = MakeMatrix(pRadiusDict, dMaxRadius, pWellFC, sWhere,
pPntDict)

If Not m_TC.Continue Then Exit Function

\' make a collection of trapezoids based on neighbor relationships

\' defined in the matrix. Trapezoids contained entirely within the old

\' unsmoothed outline are tossed.

Dim cTrapezoids As Collection

Set cTrapezoids = MakeTrapezoids2(pPntDict, pOIDDict, pRadiusDict, \_

pWellFC, sWhere, pBufferFeat.ShapeCopy)

If Not m_TC.Continue Then Exit Function

Debug.Print cTrapezoids.Count & \" trapezoids created\"

\' union the collection of trapezoids with the original field outline

\' and return that geometry

Set MakeSubFieldOutline = UnionTrapezoids(cTrapezoids, \_

pBufferFeat.ShapeCopy)

Exit Function

EH:

MsgBox \"error in MakeFieldOutlines: \" & Err.Description & \", \" &
Err.Number

End Function

Sub GetRadii(pFC As IFeatureClass, ByVal sWhere As String, \_

ByRef pDict As Scripting.Dictionary, \_

ByRef dMaxRadius As Double)

\'

\' fill up pDict with radii from the featureclass\'s \"buffer\" field

\'

\' (return the max radius also.)

On Error GoTo EH

dMaxRadius = -1

Set pDict = New Scripting.Dictionary

Dim pQF As IQueryFilter

Set pQF = New QueryFilter

pQF.WhereClause = sWhere

Dim lFld As Long

lFld = pFC.FindField(\"buffer\")

If lFld = -1 Then

MsgBox \"GetMaxRadius cannot find buffer column\"

Exit Sub

End If

Dim pFCur As IFeatureCursor

Set pFCur = pFC.Search(pQF, True)

Dim pFeat As IFeature

Set pFeat = pFCur.NextFeature

Do Until pFeat Is Nothing

Dim dRadius As Double

\' convert feet to meters

dRadius = pFeat.Value(lFld) \* 0.3048

If dRadius \> dMaxRadius Then

dMaxRadius = dRadius

End If

If Not m_TC.Continue Then Exit Sub

pDict.Add pFeat.OID, dRadius

Set pFeat = pFCur.NextFeature

Loop

Exit Sub

EH:

MsgBox \"error in GetRadii: \" & Err.Description & \", \" & Err.Number

End Sub

Function MakeMatrix(ByVal pRadiusDict As Scripting.Dictionary, \_

ByVal dMaxRadius As Double, \_

ByVal pFC As IFeatureClass, \_

ByVal sWhere As String, \_

ByRef pPntDict As Scripting.Dictionary) As Scripting.Dictionary

\'

\' populate pOIDDict keys of OIDs and list of

\' features that are within dMaxRadius

\' also populate pPntDict with the points

Dim pOIDDict As Scripting.Dictionary

Set pOIDDict = New Scripting.Dictionary

\'

\' e.g. pOIDDict.Item(22) might contain \"32,34,25\" which

\' means that features whose oids are 32,34 and 25 are within

\' 2.5x avgradius of feature whose OID is 22.

\'

\' pRadiusDict is already populated and used as a lookup

\' table for finding radii.

\'

On Error GoTo EH

Dim pQF As IQueryFilter

Set pQF = New QueryFilter

pQF.WhereClause = sWhere

Dim lFld As Long

lFld = pFC.FindField(\"buffer\")

If lFld = -1 Then

MsgBox \"makematrix cannot find buffer column\"

Exit Function

End If

Dim pFeatIndex As IFeatureIndex2

Set pFeatIndex = New FeatureIndex

Set pFeatIndex.FeatureClass = pFC

Set pFeatIndex.FeatureCursor = pFC.Search(pQF, False)

Application.StatusBar.Message(0) = \"indexing \" & pFC.AliasName

pFeatIndex.Index Nothing, GetExtent(pFC)

Dim pIQ As IIndexQuery2

Set pIQ = pFeatIndex

Application.StatusBar.ShowProgressBar \"building matrix\", 0, \_

pRadiusDict.Count, 1, True

Application.StatusBar.ProgressBar.Position = 0

Set pOIDDict = New Scripting.Dictionary

Set pPntDict = LoadDictionary(pFC, sWhere)

Dim pFCur As IFeatureCursor

Set pFCur = pFC.Search(pQF, False)

Dim pFeat As IFeature

Set pFeat = pFCur.NextFeature

Do Until pFeat Is Nothing

Dim pTopoOp As ITopologicalOperator

Set pTopoOp = pFeat.Shape

Dim pPoint1 As IPoint

Set pPoint1 = pFeat.Shape

Dim dRad1 As Double

dRad1 = pRadiusDict.Item(pFeat.OID)

Dim c As Collection

Set c = New Collection

Dim vOIDs As Variant

pIQ.IntersectedFeatures pTopoOp.Buffer(dMaxRadius \* 3#), vOIDs

Dim vOID As Variant

For Each vOID In vOIDs

Dim pPoint2 As IPoint

Set pPoint2 = pPntDict.Item(vOID)

Dim dRad2 As Double

dRad2 = pRadiusDict.Item(vOID)

Dim dAvgRadius As Double

dAvgRadius = (dRad1 + dRad2) / 2#

Dim dist As Double

dist = GetDist(pPoint1, pPoint2)

If dist \<= (3# \* dAvgRadius) And vOID \<\> pFeat.OID Then

c.Add vOID

End If

Next vOID

pOIDDict.Add pFeat.OID, c

\' Debug.Print pFeat.OID & \" has \" & c.Count & \" neighbors\"

Application.StatusBar.ProgressBar.Step

If Not m_TC.Continue Then Exit Function

Set pFeat = pFCur.NextFeature

Loop

Application.StatusBar.HideProgressBar

Debug.Print pPntDict.Count & \" points in memory\"

Set MakeMatrix = pOIDDict

Exit Function

EH:

MsgBox \"error in MakeMatrix: \" & Err.Description & \", \" & Err.Number

End Function

Function LoadDictionary(pFC As IFeatureClass, sWhere As String) As
Scripting.Dictionary

Dim pPntDict As Scripting.Dictionary

Set pPntDict = New Scripting.Dictionary

Dim pQF As IQueryFilter

Set pQF = New QueryFilter

pQF.WhereClause = sWhere

Dim pFCur As IFeatureCursor

Set pFCur = pFC.Search(pQF, False)

Dim pFeat As IFeature

Set pFeat = pFCur.NextFeature

Do Until pFeat Is Nothing

pPntDict.Add pFeat.OID, pFeat.ShapeCopy

Set pFeat = pFCur.NextFeature

Loop

Set LoadDictionary = pPntDict

End Function

Function MakeTrapezoids2(ByVal pPntDict As Scripting.Dictionary, \_

ByVal pOIDDict As Scripting.Dictionary, \_

ByVal pRadiusDict As Scripting.Dictionary, \_

ByVal pFC As IFeatureClass, \_

ByVal sWhere As String, \_

ByVal pPolygon As IPolygon) As Collection

\'

\' returns a collection of trapezoids that fill in the gaps between

\' closely neighboring wells.

\'

\'

\' create trapezoids based on the matrix (pOIDDict)

\' which describes what OIDs are near other OIDs

\'

Simplify pPolygon

\' pFC is the wells featureclass.

\' Create a spatial index that allows us to more quickly determine

\' if a trapezoid falls completely within pPolygon, in which case

\' we can omit it.

Dim pSpatialIndex As ISpatialIndex

Set pSpatialIndex = pPolygon

pSpatialIndex.AllowIndexing = True

pSpatialIndex.Invalidate

Dim pRelOp As IRelationalOperator

Set pRelOp = pSpatialIndex

On Error GoTo EH

Set MakeTrapezoids2 = New Collection

Application.StatusBar.ShowProgressBar \"making trapezoids\", 0, \_

pPntDict.Count, 1, True

Application.StatusBar.ProgressBar.Position = 0

Dim pQF As IQueryFilter

Set pQF = New QueryFilter

pQF.WhereClause = sWhere

Dim pFCur As IFeatureCursor, pFeat As IFeature

Set pFCur = pFC.Search(pQF, False)

\' loop through each well point

Set pFeat = pFCur.NextFeature

Do Until pFeat Is Nothing

Dim pPoint1 As IPoint

Set pPoint1 = pPntDict.Item(pFeat.OID)

Dim dRadius1 As Double

dRadius1 = pRadiusDict.Item(pFeat.OID)

Dim cNeighbors As Collection

Set cNeighbors = pOIDDict.Item(pFeat.OID)

\' Debug.Print pFeat.OID & \" has \" & cNeighbors.Count & \" neighbors\"

If pFeat.OID = 80 Then

Debug.Print \"ok\"

End If

Dim l As Long

For l = 1 To cNeighbors.Count

Dim lOID As Long

lOID = cNeighbors.Item(l)

If lOID \> pFeat.OID Then

Dim pTrapezoid As IPolygon

Set pTrapezoid = \_

MakeTrapezoid3(pPntDict.Item(pFeat.OID), pRadiusDict.Item(pFeat.OID), \_

pPntDict.Item(lOID), pRadiusDict.Item(lOID))

If Not pTrapezoid Is Nothing Then

If Not pRelOp.Contains(pTrapezoid.Envelope) Then

If Not pRelOp.Contains(pTrapezoid) Then

Application.StatusBar.ProgressBar.Message = \_

\"making trapezoids (\" & MakeTrapezoids2.Count & \")\"

\' add the trapezoid to the collection for later use.

MakeTrapezoids2.Add pTrapezoid

End If

Else

\'Debug.Print \"skipping \"

\' the trapezoid falls completely within the original

\' field outline, so there\'s no need to keep it.

End If

End If

Else

\' its a symmetric matrix, so don\'t worry about

\' this half (otherwise duplicates would be created

End If

If Not m_TC.Continue Then Exit Function

Next l

Application.StatusBar.ProgressBar.Step

Set pFeat = pFCur.NextFeature

Loop

Application.StatusBar.HideProgressBar

Exit Function

EH:

MsgBox \"error in MakeTrapezoids2: \" & Err.Description & \", \" &
Err.Number

End Function

Function MakeTrapezoid3(pPoint1 As IPoint, dRadius1 As Double, \_

pPoint2 As IPoint, dRadius2 As Double) As IPolygon

\'

\' Makes a trapezoid that connects two wells, based on how close

\' they need to be.

\'

Dim pOutTrapezoid As IPolygon

Set pOutTrapezoid = Nothing

Dim dAvgRadius As Double

dAvgRadius = (dRadius1 + dRadius2) / 2#

Dim dDist As Double

dDist = GetDist(pPoint1, pPoint2)

\' Jack says:

\' \"How about using the full diameter trapezoid for all

\' pairs with \< =2 distance and the ½ for all

\' from \>2 to 2.5. So those buffers that overlap or touch

\' would get a tangent from the full diameter and those

\' that are separated (but \<= 2.5) would get a \'neck\'

\' based on a smaller circle (radii).\"

If dDist \< (dAvgRadius \* 2.5) And dDist \> 0# Then

Dim pTrapezoid As IPolygon

Dim dRadiusFactor As Double

If dDist \> (2# \* dAvgRadius) Then

dRadiusFactor = 0.5

Else

dRadiusFactor = 1#

End If

\' make a trapezoid based on the two points and the radii

Set pTrapezoid = MakeTrapezoid(pPoint1, dRadiusFactor \* dRadius1, \_

pPoint2, dRadiusFactor \* dRadius2)

Simplify pTrapezoid

Set pOutTrapezoid = pTrapezoid

Else

\' the two points are close, but not close enough

\' no trapezoid created

End If

Set MakeTrapezoid3 = pOutTrapezoid

End Function

\'Function UnionTrapezoids(cTrapezoids As Collection, \_

\' pPolygon As IPolygon) As IPolygon

\' On Error GoTo EH

\' Dim pWSF As IWorkspaceFactory

\' Set pWSF = New FileGDBWorkspaceFactory

\'

\' Dim pFWS As IFeatureWorkspace

\' Set pFWS =
pWSF.OpenFromFile(\"C:\\Projects\\EIA\\Kirk_Smooth3\\Scratch.gdb\", 0)

\'

\' Dim pDS As IDataset

\' On Error Resume Next

\' Set pDS = pFWS.OpenFeatureClass(\"Trapezoids\")

\' If Err.Number = 0 Then

\' pDS.Delete

\' End If

\' Err.Clear

\' On Error GoTo EH

\'

\' Dim pFlds As IFieldsEdit

\' Set pFlds = New esrigeodatabase.Fields

\' pFlds.AddField MakeField(\"ObjectID\", esriFieldTypeOID)

\' pFlds.AddField
m_outFC.Fields.Field(m_outFC.FindField(m_outFC.ShapeFieldName))

\'

\' Dim pFC As IFeatureClass

\' Set pFC = pFWS.CreateFeatureClass(\"Trapezoids\", pFlds,
m_outFC.CLSID, Nothing, esriFTSimple, \"Shape\", \"\")

\'

\' Dim pFCur As IFeatureCursor

\' Set pFCur = pFC.Insert(False)

\'

\' Dim l As Long

\' For l = 1 To cTrapezoids.Count

\' Dim pFBuff As IFeatureBuffer

\' Set pFBuff = pFC.CreateFeatureBuffer

\' Set pFBuff.Shape = cTrapezoids.Item(l)

\' pFCur.InsertFeature pFBuff

\' Next l

\' pFCur.Flush

\' Exit Function

\'EH:

\' MsgBox \"error in UnionTrapezoids: \" & Err.Description & \", \" &
Err.Number

\'End Function

Function UnionTrapezoids(cTrapezoids As Collection, \_

pPolygon As IPolygon) As IPolygon

On Error GoTo EH

Application.StatusBar.ShowProgressBar \"unioning trapezoids\", 0, \_

cTrapezoids.Count, 1, True

Application.StatusBar.ProgressBar.Position = 0

\' ITopologicalOperator.Union was found to raise an error

\' occasionally, so the code below relies on a workaround.

\'

\' The workaround involves keeping a list of bad trapezoids and

\' also a backup polygon. When a trapezoid causes an error, it

\' is added to the list and processing starts over at the last

\' savepoint.

\'

Simplify pPolygon

Dim pBackup As IPolygon

Set pBackup = Clone(pPolygon)

Dim pBadOnes As New Scripting.Dictionary

Dim pTopoOp As ITopologicalOperator

Set pTopoOp = Clone(pPolygon)

Dim l As Long

l = 1

Dim lBackupID As Long

lBackupID = 1

Do While l \< cTrapezoids.Count + 1

Dim pTrapezoid As IPolygon

Set pTrapezoid = cTrapezoids.Item(l)

If pBadOnes.Exists(l) Then

Debug.Print \"moving bad trapezoid\" & l & \_

\" attempt \" & pBadOnes.Item(l)

Dim pTrans2D As ITransform2D

Set pTrans2D = pTrapezoid

\' move 10 centimeters

pTrans2D.Move 0.1, 0.1

End If

On Error Resume Next

Set pTopoOp = pTopoOp.Union(pTrapezoid)

If Err.Number = 0 Then

l = l + 1

\' every 100 steps make a backup

If l Mod 100 = 0 Then

Set pBackup = Clone(pTopoOp)

lBackupID = l

Debug.Print l

End If

Application.StatusBar.ProgressBar.Step

Else

Application.StatusBar.ProgressBar.Message = \_

\"errorcount: \" & pBadOnes.Count

If Not pBadOnes.Exists(l) Then

pBadOnes.Add l, 1

Else

\' a repeat offender !

Debug.Print \"repeat offender encountered\"

pBadOnes.Item(l) = pBadOnes.Item(l) + 1

End If

l = lBackupID

Set pTopoOp = pBackup

Application.StatusBar.ProgressBar.Position = l

End If

If Not m_TC.Continue Then Exit Function

On Error GoTo EH

Loop

Application.StatusBar.HideProgressBar

Set UnionTrapezoids = pTopoOp

Exit Function

EH:

MsgBox \"error in UnionTrapezoids: \" & Err.Description & \", \" &
Err.Number

End Function

Function GetExtent(pGDS As IGeoDataset) As IEnvelope

Set GetExtent = pGDS.Extent

End Function

Function GetDist(pProxOp As IProximityOperator, pGeom As IGeometry) As
Double

Dim pPCS As IProjectedCoordinateSystem

Set pPCS = pGeom.SpatialReference

Dim d As Double

d = pProxOp.ReturnDistance(pGeom) \* pPCS.CoordinateUnit.MetersPerUnit

GetDist = d

End Function

Function MakeTrapezoid(pPoint1 As IPoint, dDist1 As Double, \_

pPoint2 As IPoint, dDist2 As Double) As IPointCollection

On Error GoTo EH

\' create a new polygon based on points and distances

Set MakeTrapezoid = New Polygon

Dim pLine As ILine

Set pLine = New esriGeometry.Line

pLine.PutCoords pPoint1, pPoint2

Dim pCPoint As IConstructPoint, dDist As Double

Set pCPoint = New Point

pCPoint.ConstructOffset pLine, esriNoExtension, 0#, True, dDist1

MakeTrapezoid.AddPoint pCPoint

pCPoint.ConstructOffset pLine, esriNoExtension, 1#, True, dDist2

MakeTrapezoid.AddPoint pCPoint

pCPoint.ConstructOffset pLine, esriNoExtension, 1#, True, -dDist2

MakeTrapezoid.AddPoint pCPoint

pCPoint.ConstructOffset pLine, esriNoExtension, 0#, True, -dDist1

MakeTrapezoid.AddPoint pCPoint

Exit Function

EH:

MsgBox \"error in MakeTrapezoid: \" & Err.Description & \", \" &
Err.Number

End Function

Sub Simplify(pTopoOp2 As ITopologicalOperator2)

pTopoOp2.IsKnownSimple = False

pTopoOp2.Simplify

End Sub

Function Clone(pClone As IClone) As IClone

Set Clone = pClone.Clone

End Function

Function GetArea(pArea As IArea) As Double

GetArea = pArea.Area

End Function

Sub CopyColumns(pFBuff As IFeatureBuffer, pInFeat As IFeature)

\' copy all field values from pInFeat to pFBuff based

\' on field names. Skip the geometry and non-editable fields

Dim lOutFld As Long

For lOutFld = 0 To pFBuff.Fields.FieldCount - 1

If pFBuff.Fields.Field(lOutFld).Editable Then

If pFBuff.Fields.Field(lOutFld).Type \_

\<\> esriFieldTypeGeometry Then

Dim lInFld As Long

lInFld = pInFeat.Fields. \_

FindField(pFBuff.Fields.Field(lOutFld).Name)

If lInFld \> -1 Then

pFBuff.Value(lOutFld) = pInFeat.Value(lInFld)

End If

End If

End If

Next lOutFld

End Sub

Function MakeField(sName As String, lType As esriFieldType, \_

Optional lLen As Long) As esrigeodatabase.IField

\' make a new field

Dim pFldEdit As IFieldEdit

Set pFldEdit = New esrigeodatabase.Field

pFldEdit.Name = sName

pFldEdit.Type = lType

If lLen \> 0 Then

pFldEdit.Length = lLen

End If

Set MakeField = pFldEdit

End Function

Function CreateFC(pInFC As IFeatureClass, sName As String) As
IFeatureClass

\'

\' create a featureclass identical to pInFC and in the

\' same workspace as pInFC, but named sName

\'

Dim pFWS As IFeatureWorkspace

Set pFWS = GetWS(pInFC)

Dim pDS As IDataset

On Error Resume Next

Set pDS = pFWS.OpenFeatureClass(sName)

If Err.Number = 0 Then

Debug.Print \"deleting \" & pDS.Name

pDS.Delete

Else

Err.Clear

End If

On Error GoTo EH

Dim pFldsEdit As IFieldsEdit

Set pFldsEdit = Clone(pInFC.Fields)

If pInFC.FindField(\"PctChange\") = -1 Then

pFldsEdit.AddField MakeField(\"PctChange\", esriFieldTypeDouble)

End If

Set CreateFC = pFWS.CreateFeatureClass(sName, pFldsEdit, \_

pInFC.CLSID, Nothing, \_

esriFTSimple, \_

pInFC.ShapeFieldName, \"\")

Exit Function

EH:

MsgBox \"error in CreateFC: \" & Err.Description & \", \" & Err.Number

End Function

Function GetWS(pDS As IDataset) As IWorkspace

Set GetWS = pDS.Workspace

End Function