评估网格化结果

原始讨论：http://goo.gl/VlzzNT

评估网格化结果的一般方法

对原始数据网格化之后，生成等值线图，叠加原始数据点位图并标注数据，是最直接也是最常用的方法。

如果要写评估报告，就还要更进一步，用上一些数值方法，获得一些评估参数。

比如分别计算原始数据和残差的最小值、最大值、平均值、标准误差、标准离差等，然后比对，就可评估网格化结果与原始数据的吻合情况。

再如，用一个叫拟合优度（Goodness of Fit）的参数，可以非常直接地指示网格化的精度。如果拟合优度为 0，表示网格化模型与原始数据完全无关；如果拟合优度为 1，表示网格化模型与原始数据每个点都精确匹配。所以，拟合优度越高，网格化结果越理想。

最后，如果要对网格化效果不佳的区域进行数据再处理，例如增加数据点。可以绘制残差分布等值线图，叠加数据点位图层，即可直观、全面地查看网格化情况。

实现网格结果评估自动化

'=====================================================
'Numerically Verify Grid Results
'Golden Software Surfer 13.x + Scripter
'Holz
'2016-5-4
'=====================================================
Option Explicit
Sub Main
    Dim SurferApp As Object
    Dim Plot As Object
    Dim wks As Object
    Dim stats As Object
    Dim xCol As Integer, yCol As Integer, zCol As Integer
    Dim rCol As Integer
    Dim datFile As String, grdFile As String

    xCol = 1
    yCol = 2
    zCol = 3
    rCol = 4

    Debug.Clear
    '已知当前目录下有原始数据文件
    datFile = CurDir() & "\demogrid.dat"
    '如果没有，则弹出对话框选择原始数据文件
    If Dir(datFile) = "" Then
        datFile = GetFilePath(,"csv|dat|txt",,"选择原始数据",4)
    End If
    '如果仍然未选择任何文件，程序结束
    If datFile = "" Then End
    Debug.Print "======= " & Date & " " & Time & " ======="
    Debug.Print "原始数据： " & datFile

    '已知当前目录下有原始数据文件对应的网格文件
    grdFile = CurDir() & "\demogrid.grd"
    '如果没有，则弹出对话框选择原始数据文件
    If Dir(grdFile) = "" Then
        grdFile = GetFilePath(,"grd",,"选择原始数据对应的网格文件",4)
    End If
    '如果仍然未选择任何文件，程序结束
    If grdFile = "" Then End
    Debug.Print "对应网格： " & grdFile

    Set SurferApp = CreateObject("Surfer.Application")
    Set wks = SurferApp.Documents.open(datFile)
    Dim r1 As Integer, r2 As Integer
    r1 = 2
    r2 = wks.UsedRange.RowCount
    rCol = wks.UsedRange.ColumnCount + 1
    Debug.Print "末尾行号： " & r2
    Debug.Print "残差列号： " & rCol
    wks.Close(srfSaveChangesNo)

    '计算残差
    Set wks = SurferApp.GridResiduals2(grdFile, datFile, _
        xCol, yCol, zCol, rCol)
    wks.Cells(1,rCol).Value = "GridResiduals"
    '用统计方法计算原始数据的最小、最大、均值、标准差、离差等
    Set stats = wks.Columns(zCol).Statistics(True, True, wksStatsAll)
    Dim zMin As Double, zMax As Double, zMean As Double
    Dim zSE As Double, zSD As Double
    zMin = stats.Minimum
    zMax = stats.Maximum
    zMean = stats.Mean
    zSE = stats.StandardError
    zSD = stats.StandardDeviation
    '用统计方法计算残差的最小、最大、均值、标准差、离差等
    Set stats = wks.Columns(rCol).Statistics(True, True, wksStatsAll)
    Dim rMin As Double, rMax As Double, rMean As Double
    Dim rSE As Double, rSD As Double
    rMin = stats.Minimum
    rMax = stats.Maximum
    rMean = stats.Mean
    rSE = stats.StandardError
    rSD = stats.StandardDeviation
    '通过原始数据与残差的统计值比对，可评估网格化的精度。
    '比如，残差最大值与原始数据范围之间的比值越小，说明
    '变异小，网格化方法和参数选择得当，网格化结果忠于数据。
    '残差平均值越接近于0，说明网格化结果越接近于原始数据。
    '标准离差的大小指示了部分有偏差的区域，其网格化结果
    '偏离数据点的情况。
    Debug.Print "------------------------------------"
    Debug.Print "统计项目 " & vbTab & "原始数据" & vbTab & "残差"
    Debug.Print "  最小值 " & vbTab & Round(zMin,4) & vbTab & vbTab & Round(rMin,4)
    Debug.Print "  最大值 " & vbTab & Round(zMax,4) & vbTab & vbTab & Round(rMax,4)
    Debug.Print "  平均值 " & vbTab & Round(zMean,4) & vbTab & vbTab & Round(rMean,4)
    Debug.Print "  标准差 " & vbTab & Round(zSE,4) & vbTab & vbTab & Round(rSE,4)
    Debug.Print "标准离差 " & vbTab & Round(zSD,4) & vbTab & vbTab & Round(rSD,4)
    Debug.Print "------------------------------------"

    Dim r2Col As Integer
    r2Col = rCol + 1
    Dim eq As String
    eq = Chr(64+rCol)
    eq = Chr(64+r2Col) & "=" & eq & "*" & eq
    Debug.Print "残差平方公式： " & eq
    wks.Transform3(wksTransformWithCol, r1, r2, eq)
    Dim rSum As Double
    rSum = wks.Columns(r2Col).Statistics(wksStatsSum).Sum
    Debug.Print "  残差平方和： " & Round(rSum,4)
    eq = "(" & zMean & "-" & Chr(64+zCol) & ")"
    eq = Chr(64+r2Col) & "=" & eq & "*" & eq
    Debug.Print "Z 值均差平方公式： " & eq
    wks.transform2(wksTransformWithCol, r1, r2, eq)
    Dim zSum As Double
    zSum = wks.Columns(r2Col).Statistics(wksStatsSum).Sum
    Debug.Print "  Z 值均差平方和： " & Round(zSum,4)
    Dim FoG As Double
    FoG = 1.0 - rSum / zSum
    '拟合优度（Goodness of Fit）非常直接地指示了网格化的精度。
    '如果拟合优度为 0，表示网格化模型与原始数据完全无关；
    '如果拟合优度为 1，表示网格化模型与原始数据每个点都精确匹配。
    '所以，拟合优度越高，网格化结果越理想。
    Debug.Print "拟合优度： " & Format(FoG,"0.0000")
    wks.Columns(r2Col).Delete(wksDeleteColumns)
    wks.Close(srfSaveChangesYes)

    '用克里格法网格化残差列，得到残差网格
    grdFile = Left(datFile,InStrRev(datFile,".")-1) & "_r.grd"
    Dim bstatus As Boolean
    bstatus = SurferApp.GridData3(DataFile:=datFile, _
        OutGrid:=grdFile, ShowReport:=False, Algorithm:=2, _
        xCol:=xCol, yCol:=yCol, zCol:=rCol)
    Debug.Print "  网格化数据： " & If(bstatus,"","不") & "成功。"
    Set wks = SurferApp.Documents.open(datFile)
    wks.Columns(rCol).Delete(wksDeleteColumns)
    wks.Close(srfSaveChangesYes)

    '新开一个绘图文档
    Set Plot = SurferApp.Documents.Add(srfDocPlot)
    Dim MapFrame As Object
    Dim MapLayer As Object
    '用刚生成的残差网格生成等值线图
    Set MapFrame = Plot.Shapes.AddContourMap(grdFile)
    '设置地图的大小、比例
    MapFrame.xLength = 27
    MapFrame.yMapPerPU = MapFrame.xMapPerPU
    '设置地图的坐标参数
    Dim Axis As Object
    For Each Axis In MapFrame.Axes
        Axis.MajorTickType = srfTickNone
        Axis.MinorTickType = srfTickNone
        Axis.ShowLabels = False
        '因为标题用了中文，必须使用中文字体
        Axis.TitleFont.Face = "宋体"
        Axis.TitleFont.ForeColorRGBA.Color = srfColorBlue
        If Axis.AxisType = srfATTop Then
            Axis.Title = "残差分布等值线图"
            Axis.TitleFont.Size = 40
        End If
        If Axis.AxisType = srfATBottom Then
            '可以在字符串中使用数学文本指令
            Axis.Title = "拟合优度 R\up50 {\fs50 2}\dn50  = " & Format(FoG,"0.0000")
            Axis.TitleFont.Size = 20
        End If
    Next

    '设置等值线图层
    Set MapLayer = MapFrame.Overlays(1)
    With MapLayer
        'Surfer 13 的等级设置有许多联动参数，
        '所以先将等级方法设置为简单。
        .LevelMethod = SrfConLevelMethodSimple
        .SmoothContours = srfConSmoothMed
        .MajorLine.Style = "Invisible"
        .MinorLine.Style = "Invisible"
        .ShowMajorLabels = False
        .ShowMinorLabels = False
        '要设置等值线填充，第一步开启填充
        .FillContours = True
        Dim pos(2) As Double
        pos(0) = 0.0
        pos(1) = .FillForegroundColorMap.DatToPos(0)
        pos(2) = 1.0
        Dim colors(2) As Long
        colors(0) = srfColorRed
        colors(1) = srfColorWhite
        colors(2) = srfColorRed
        '第二步载入或设置色谱。
        'Surfer 13 帮助文件里居然没有 SetNodes
        .FillForegroundColorMap.SetNodes(Positions:=pos, colors:=colors)
        '第三步应用填充设置，不按部就班就得不到结果。
        .ApplyFillToLevels(1, 1, 0)
        .ShowColorScale = True
        .ColorScale.FirstLabel = 2
    End With

    '添加一个数据点位图层，配合等值线图，就可以直观地
    '看到哪些区域与原始数据最吻合，哪些区域偏差大。
    '对照此图就可以决定哪里需要加密数据，增加数据。
    Set MapLayer = Plot.Shapes.AddPostLayer(MapFrame, _
        datFile, xCol, yCol)

    SurferApp.ActiveWindow.Zoom(srfZoomFitToWindow)
    SurferApp.Visible = True
    Debug.Print "======= " & Date & " " & Time & " ======="
End Sub