After surveying and mapping a large area, it can sometimes be necessary to calculate statistics such as average elevation or peak contamination within multiple smaller areas of the map.  The script below shows one way this can be done via automation.  In this script

• an SHP file of township & range polygons is exported to individual BLN files
• the grid file (contours) is assigned NoData outside each individual BLN and converted to DAT.
• the average Z value of each DAT is calculated and saved to a new worksheet.

The end result is a new worksheet containing the average Z value (contour) within each polygon (township or range).

To run this script:

1. Copy the script below or download the attached Average Z in Polygon.bas.
2. In a Windows Explorer window, navigate to C:\Program Files\Golden Software\Surfer.
3. Double click on Scripter.exe to launch Scripter.
4. Press Ctrl+A to select all of the existing lines then press Delete.
5. If you copied this script, press Ctrl+V to paste it into Scripter. If you downloaded it, click File | Open, select the BAS file from your downloads directory, and click Open.
6. Update the values on lines 9, 13, 17, and 70 as directed
7. Click Script | Run to run the script.
    

*********

Sub Main
    'Initializes and opens Surfer 
    Dim SurferApp, Plot, MapFrame, BaseMap, ContourMap, Wks, WksRange, Wks2, WksRange2, WksRange3 As Object
    Set SurferApp = CreateObject("Surfer.Application")
    SurferApp.Visible = True
    Set Plot = SurferApp.Documents.Add(srfDocPlot)

    'Sets the path 
    Path$="C:\users\gsadmin\desktop\" '*****Here type in the file path containing your GRD and SHP files******** 
    OutPath$=Path$+"BLN\"

    'Create base map of Township/Range polygons 
    Set MapFrame = Plot.Shapes.AddBaseMap(ImportFileName:=Path$+"boundarypoly.shp") '*****Here type in name of your SHP file******** 
    Set BaseMap = MapFrame.Overlays(1)

    'Exports all of the polygons to a BLN 
    numpoly = 19 '*****Here type in the number of polygons that are in your township/range SHP file******** 
    Dim Header (19) As Integer
    MapFrame.Axes(1).Visible=False
    MapFrame.Axes(2).Visible=False
    MapFrame.Axes(3).Visible=False
    MapFrame.Axes(4).Visible=False
    Plot.Export(FileName:=Path$+"township_range.bln", SelectionOnly:=False)

    'Opens the BLN and splits the polygons out into separate BLN files 
    Set Wks = SurferApp.Documents.Open(FileName:=Path$+"township_range.bln")
    Set WksRange = Wks.Columns("C")

        'Search through each row and save the ones without z values as header rows 
        hdrct = 1
        numrows = WksRange.RowCount
        For i=1 To numrows
            If Wks.Cells((i),3)="" Then 
                Header(hdrct) = i
                hdrct = hdrct+1
            End If
        Next
        hdrct = hdrct-1

        'Select the values from each header row to the row before the next header row and c/p those cells into a new wks. Save the wks. 
        For i=1 To hdrct-1
            Set WksRange2 = Wks.Cells(Row:=Header(i), Col:=1, LastRow:=Header(i+1)-1, LastCol:=3)
            WksRange2.Copy
            Set Wks2 = SurferApp.Documents.Add(DocType:=srfDocWks)
            Set WksRange3 = Wks2.Cells(Row:=1, Col:=1, LastRow:=Header(i+1)-Header(i), LastCol:=3)
            WksRange3.Paste
            Wks2.SaveAs(FileName:=OutPath$+i+".bln", FileFormat:=srfSaveFormatBln)
            Wks2.Close(SaveChanges:=srfSaveChangesNo)
        Next

        'Repeats what was done above for the last polygon (which doesn't have a next header to go to) 
        Set WksRange2 = Wks.Cells(Row:=Header(hdrct), Col:=1, LastRow:=numrows, LastCol:=3)
        WksRange2.Copy
        Set Wks2 = SurferApp.Documents.Add(DocType:=srfDocWks)
        Set WksRange3 = Wks2.Cells(Row:=1, Col:=1, LastRow:=numrows-Header(hdrct), LastCol:=3)
        WksRange3.Paste
        Wks2.SaveAs(FileName:=OutPath$+i+".bln", FileFormat:=srfSaveFormatBln)
        Wks2.Close(SaveChanges:=srfSaveChangesNo)

    Wks.Close(SaveChanges:=srfSaveChangesNo)

    'Creates a file to contain the results 
    Set Wks = SurferApp.Documents.Add(DocType:=srfDocWks)

    'Loops through the BLN files 
    filecount=1
    file = Dir (OutPath$+"*.bln")
    While file <> ""
        'Blank outside of BLN file 
        SurferApp.GridBlank(InGrid:=Path$+"Colorado.grd", BlankFile:=OutPath$+file, OutGrid:=OutPath$+file+".dat", OutFmt:=srfGridFmtXYZ) '*****Here type in the name of your GRD file******** 

        'Open DAT in worksheet 
        Set Wks2 = SurferApp.Documents.Open(FileName:=OutPath$+file+".dat")
        Set WksRange2 = Wks2.Columns("C")

        Set WksRange3 = Wks2.Columns(Col1:=1, Col2:=3)
        WksRange3.Sort(Col1:=3, Order1:=wksSortAscending)

        'Delete all lines where z is blanking value 
        For i=1 To WksRange2.RowCount
            If WksRange2.Cells((i))=1.70141E+38 Then 
                WksRange2.Cells(Row:=i, Col:=1, LastRow:=WksRange2.RowCount, LastCol:=3).Delete(Direction:=wksDeleteRows)
                Exit For
            End If
        Next

        'Calculate statistics (average) on z column 
        Set Stats = WksRange2.Statistics(Flags:=wksStatsMean)

        'Write Township/Range name and z average value to a file (or an array to print out) 
        Set WksRange = Wks.Columns(Col1:=1, Col2:=2)
        WksRange.Cells((filecount),1) = filecount
        WksRange.Cells((filecount),2) = Stats.Mean
        Wks.SaveAs(FileName:=OutPath$+"Averages.xls", FileFormat:=srfSaveFormatExcel)
        filecount=filecount+1

        'Get next file file = Dir() 
    Wend

End Sub

Updated November 2021