Grid Analysis using Visual Basic for Applications Exercise

Prepared by David Tarboton, Utah State University.

Purpose

The purpose of this exercise is to learn how to use Visual Basic for Application (VBA) to Extend ArcGIS.

Exercise

It is not always possible to do what is needed using the functionality provided by ArcGIS.  Therefore the capability to extend ArcGIS through programming or scripting capability has been developed using Visual Basic for Applications (VBA).  This is similar to the concept of Macro's in Excel.  Here you will implement a simple VBA Macro to calculate the flow path distance from each point in the watershed to the outlet.  [This is something that can actually be achieved using the flowpath command in ArcInfo command line GRID functionality, and perhaps through the spatial analyst object library but programming it is a good learning exercise.]

 

This exercise will illustrate how to use ArcObjects from a VBA macro that I have written for you.  The macro and demo data are in the file vbaex.zip.

 

From within ArcMAP select Tools/Macros/Macros...

Enter a Macro name: dist and click Create.

Cut and paste the following VBA code into the Visual Basic Editor that opens.  Be careful to paste code above, inbetween and below the stub code

Sub dist()

 

End Sub

that appears.  In Visual Basic, comments are indicated starting with a ', and I have used these to try explain what is going on.  Some comments are lengthy and may wrap around on your display.  These appear as red in the VB editor and after pasting you need to unwrap them, or else the script will not run. You could also copy equivalent code from vba.mxd in the vbaex.zip example file.

' BEGINNING OF CODE

Option Explicit   ' This is a command saying that all variables need to be explicitly defined.  It is good programming policy

'  Global variables.  These are declared before the first function or subroutine.  These are variables that will be used by the subroutine dist, as well as the subroutine distcalc

    Dim pPixels As Variant

    Dim dPixels

    Dim di(8) As Integer, dj(8) As Integer, dd(8) As Double

    Dim ncol As Long

    Dim nrow As Long

 

Sub dist()

' Dimensioning variables

    Dim Foldername As String

    Dim OutletShapefilename As String

    Dim pGridname As String, dGridname As String

    Dim i As Long, j As Long, k As Long

   

'*******************************************************

' Specifying inputs.  YOU WILL NEED TO CHANGE THE FOLDERNAME BELOW TO REFLECT THE LOCATION OF YOUR WORK

    Foldername = "C:\Dave\Ex7"

' YOU WILL NEED TO CHANGE THE FILENAMES BELOW TO REFLECT THE NAMES YOU ARE USING

    OutletShapefilename = "outlet"  ' THIS IS THE SHAPEFILE CONTAINING THE OUTLETS

    pGridname = "demp"             ' THIS IS THE D8 FLOW DIRECTIONS FILE CREATED BY THE COMMAND 'D8 FLOW DIRECTIONS'

    dGridname = "dist"        ' THIS IS THE NAME OF THE GRID FILE THAT WILL BE OUTPUT

'********************************************************

'  General ArcGIS objects that need to be initiated

    Dim pDoc As IMxDocument

    Set pDoc = ThisDocument

' Create the workspace factory and open shapefile feature class

   Dim pWSF As IWorkspaceFactory    ' This says that the name pWSF will refer to a WorkspaceFactory Object

   Dim pFWS As IFeatureWorkspace    ' This says that the name pFWS will refer to a Feature Workspace (in this case a folder containing feature datasets)

   Set pWSF = New ShapefileWorkspaceFactory    ' This instantiates (creates a new instance of) a WorkspaceFactory object

   Set pFWS = pWSF.OpenFromFile(Foldername, 0)   ' The function OpenFromFile that is part of (a method within) the Workspace Factory object is used to open the given folder as a Feature WorkSpace and assign the result to the named FeatureWorkspace object pFWS.

   Dim pFeClass As IFeatureClass    ' This says that the name pFeClass will refer to a Feature Class (in this case a shapefile).

  Set pFeClass = pFWS.OpenFeatureClass(OutletShapefilename) ' The function OpenFeatureClass is part of (a method within) the FeatureWorkspace object and here returns the feature class object (in this case shapefile) being worked on.

        ' OutletShapefilename is the shapefile name without suffix like "outlet"

'READ THE COORDINATES FROM THE OUTLET SHAPEFILE

' get the number of features

    Dim nFeCount As Integer

    nFeCount = pFeClass.FeatureCount(Nothing)

    Dim dbXX() As Double, dbYY() As Double  ' X and Y coordinates will be stored in double precision arrays

    ReDim dbXX(nFeCount)

    ReDim dbYY(nFeCount)

    Dim pFeature As IFeature ' A feature object variable name

    Dim pPoint As IPoint ' a point object variable name

'Loop through to get the x and y coordinates of each point associated with each feature in the feature class

    For i = 0 To nFeCount - 1

        Set pFeature = pFeClass.GetFeature(i)

        Set pPoint = pFeature.Shape 'if it is a point shapefile

        'to get the coordinate of a point

        dbXX(i) = pPoint.X

        dbYY(i) = pPoint.Y

        MsgBox "Outlet " & dbXX(i) & " " & dbYY(i)   ' This pops up a message box giving the coordinates of each outlet

    Next i

' READ THE FLOW DIRECTION GRID

' Now work with the flow direction grid

'  Raster workspace factory object

'  Dim pWSF As IWorkspaceFactory    '  This statement would generally be needed, but is not here because we are reusing the same workspace factory name used above.

    Set pWSF = New RasterWorkspaceFactory

' Raster workspace object

    Dim pRWS As IRasterWorkspace

    Set pRWS = pWSF.OpenFromFile(Foldername, 0)

' Raster dataset object

    Dim pRsDS As IRasterDataset

    Set pRsDS = pRWS.OpenRasterDataset(pGridname)

' Raster object

   Dim pRs As IRaster

   Set pRs = pRsDS.CreateDefaultRaster

' Raster band collection object.  Rasters may contain multiple bands

   Dim pRsBC As IRasterBandCollection

   Set pRsBC = pRs

' Raster band object

    Dim pRsBand As IRasterBand

    Set pRsBand = pRsBC.Item(0)    ' The first band, i.e. item number 0 in the band collection

' Raster properties object

    Dim pRsProp As IRasterProps

    Set pRsProp = pRsBand

'  Get the size and extent of flow direction grid

    Dim xll As Double

    Dim yll As Double

    Dim csize As Double

    xll = pRsProp.Extent.XMin

    yll = pRsProp.Extent.YMin

    ncol = pRsProp.Width

    nrow = pRsProp.Height

    csize = pRsProp.MeanCellSize.X

' CREATE DISTANCE TO OUTLET GRID WITH NO DATA VALUES

' Spatial reference object for new grid

    Dim pOutSR As ISpatialReference   ' Get the spatial reference of the input grid

    Set pOutSR = pRsProp.SpatialReference

'  Raster workspace object that has method to create new grid

    Dim pRWS2 As IRasterWorkspace2

    Set pRWS2 = pWSF.OpenFromFile(Foldername, 0)

' Origin object for new grid

    Dim nOrigin As IPoint

    Set nOrigin = New Point

    nOrigin.PutCoords xll, yll

' Create new grid as raster dataset object

    Dim pOutDs As IRasterDataset

    Set pOutDs = pRWS2.CreateRasterDataset(dGridname, "GRID", nOrigin, _

    ncol, nrow, csize, csize, 1, PT_FLOAT, pOutSR, True)

'  new raster indicated by "d" for distance at the beginning of the name

' new band collection

    Dim dRsBC As IRasterBandCollection

    Set dRsBC = pOutDs

' Raw Pixel objects.  These are used to access the raw data in grids

    Dim pRawPixels As IRawPixels, dRawPixels As IRawPixels

    Set pRawPixels = pRsBand

    Set dRawPixels = dRsBC.Item(0)

' new properties

    Dim dRsProp As IRasterProps

    Set dRsProp = dRawPixels

'  No data value

   Dim noDataValue As Double

   noDataValue = dRsProp.noDataValue

' Set up double point object specifying the extent of the part of the grid to work with.  In this case the entire grid.

    Dim pPnt As IPnt

    Set pPnt = New DblPnt

    pPnt.SetCoords pRsProp.Width, pRsProp.Height

' Pixel block objects to work with,within extent specified by pPnt (in this case the entire grid).

    Dim pPixelBlock As IPixelBlock, dPixelBlock As IPixelBlock3

    Set pPixelBlock = pRawPixels.CreatePixelBlock(pPnt)

    Set dPixelBlock = dRawPixels.CreatePixelBlock(pPnt)

' Origin of coordinates to be read from raw pixels into pixel block

  Dim pOrigin As IPnt

    Set pOrigin = New DblPnt

    pOrigin.SetCoords 0, 0

' Read block of flow direction pixels into PixelBlock

    pRawPixels.Read pOrigin, pPixelBlock

    dRawPixels.Read pOrigin, dPixelBlock

' Convert PixelBlock into arrays that can be referenced by rows and columns.  These are dimensioned before the sub statement so that they are global variables accessible to other modules

    pPixels = pPixelBlock.SafeArray(0)

    dPixels = dPixelBlock.PixelDataByRef(0)

' Initialize the distance pixels to no data

    For i = 0 To nrow - 1

        For j = 0 To ncol - 1

           ' col, row

            dPixels(j, i) = noDataValue

        Next j

    Next i

'PREPARE ARRAYS FOR NAVIGATING GRID

'   Set up directions

'  row offsets

    di(1) = 0

    di(2) = -1

    di(3) = -1

    di(4) = -1

    di(5) = 0

    di(6) = 1

    di(7) = 1

    di(8) = 1

'  column offsets

    dj(1) = 1

    dj(2) = 1

    dj(3) = 0

    dj(4) = -1

    dj(5) = -1

    dj(6) = -1

    dj(7) = 0

    dj(8) = 1

'  distances

    For k = 1 To 8

        dd(k) = (((csize * di(k)) ^ 2 + (csize * dj(k)) ^ 2)) ^ 0.5   ' Pythagorous Theorem

    Next k

' CONVERT EACH OUTLET COORDINATE TO ROW AND COLUMN AND START DISTANCE CALCULATION AT IT

    For k = 0 To nFeCount - 1

       i = nrow - Int((dbYY(k) - yll) / csize) - 1

       j = (dbXX(k) - xll) / csize

       dPixels(j, i) = 0

       DistCalc i, j   ' This is a recursive subroutine call

    Next k

'  Save the result

    Dim dCache

    Set dCache = dRawPixels.AcquireCache

    dRawPixels.Write pOrigin, dPixelBlock

    dRawPixels.ReturnCache dCache

   

'  Add the new raster layer to the document display

    Dim ROutLayer As IRasterLayer

    Set ROutLayer = New RasterLayer

    ROutLayer.CreateFromDataset pOutDs

    pDoc.FocusMap.AddLayer ROutLayer

    pDoc.ActiveView.Refresh

 

End Sub

 

'RECURSIVE DISTANCE CALCULATION FUNCTION

Sub DistCalc(i, j)

Dim k As Integer, inb As Long, jnb As Long

For k = 1 To 8  ' for each neighbor

    inb = i + di(k)

    jnb = j + dj(k)

    If (inb >= 0 And inb < nrow And jnb >= 0 And jnb < ncol) Then   ' guard against out of domain

        If pPixels(jnb, inb) > 0 Then   ' guard against no data

            If (pPixels(jnb, inb) - 4 = k Or pPixels(jnb, inb) + 4 = k) Then

            ' Here we have a grid cell that drains back to the grid cell we are at

                dPixels(jnb, inb) = dPixels(j, i) + dd(k)

                DistCalc inb, jnb ' Call the function for that pixel

            End If

        End If

    End If

Next k

End Sub

 

 

Once the code is entered we need some data to test it.  Switch back to ArcMap and use the Ascii to Raster tool to import the file dem.asc from vbaex.zip and save as a grid named DEM with output type float.

 

 

Use Taudem/Fill Pits to create the grid 'demfel'.  (There are no pits in this small 6x8 dataset but this is the easiest way to get the file demfel and be consistent with the naming convention.)

 

Use TauDEM/D8 Flow Directions to create the flow direction grid 'demp'. 

 

 

This will serve as the flow direction grid for the calculation of distances to the outlet.

 

Create a point shapefile with a single point within the grid cell 3rd across in the first row. 

 

 

Switch back to the Visual Basic Editor. 

Modify the code at the location about 10 lines from the top where it indicates that changes are needed to set the variables Foldername, OutletShapefilename, pGridname and dGridname to be what you are using.  Save your project.  Then run the script by clicking on the Run Sub/UserForm arrow

 

If all goes well you should have a new grid that contains distances to the outlet from each point within the watershed that drains to the outlet.  I have found that the color scheme scale is not nicely set by this procedure when the layers are added so it may be necessary to remove then add the layer.

 

 

Congratulations!  You have completed a GIS VBA program.  Now the only limit on what you can do is your creativity.

 

Use Insert/Data Frame to insert a new data frame, then remove the old one. 

 

Add the Logan River data from the Logan River Exercise to the new data frame.  Run the distance to the outlet function for the Logan River.

 

Prepare a layout that shows the distance to the outlet of the Logan River.  Include the stream network in this layout.  Report the longest distance to the outlet from within the Logan River Watershed.

These materials may be used for study, research, and education, but please credit the authors and the Utah Water Research Laboratory, Utah State University. All commercial rights reserved. Copyright 2004 Utah State University.