Determine the average contour value within each township/range polygon using Surfer automation
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:
- Copy the script below or download the attached Average Z in Polygon.bas.
- In a Windows Explorer window, navigate to C:\Program Files\Golden Software\Surfer.
- Double click on Scripter.exe to launch Scripter.
- Press Ctrl+A to select all of the existing lines then press Delete.
- 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.
- Update the values on lines 9, 13, 17, and 70 as directed
- 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
Comments
0 comments
Please sign in to leave a comment.