用 Surfer 创建热度图或密度图
本文使用 Surfer 13.x 完成。
概述
热度图（Heat Maps）或密度图（Density Maps）用来可视化特定区域内的密度区、关注区、缓冲区等。例如用热度图来显示全球的地震活动区域，虽然全球到处都有地震发生，但地震高发区（或易发区）却只是一小部分。
在 Surfer 中，所谓热度图（或密度图），只是描述了特定区域内数据出现的次数或个数，所以用 Surfer 绘制热度图其实只需要 XY 两列数据即可。
作为练手用，我从 USGS (http://earthquake.usgs.gov/earthquakes/map/) 下载了 30 天的全球地震数据，用 Surfer 把数据可视化，从而直观地查看全球范围内，都有哪些地方是地震高发区。
绘制热度图的操作步骤
选择 Grid | Data 命令启动网格化过程，在弹出的对话框中选择要用来绘制热度图的数据文件，这个文件只需要有表示 XY 位置的两列数据即可。
在网格化数据对话框中，X、Y列选择数据文件中相应的列即可；Z 列可以随意选用 X 列或 Y 列数据，因为我们只要统计点密度，所以这个 Z 值实际上没什么用处；网格化方法则要选用 Data Metrics，这个方法不是一般意义上的网格化方法，它只是一系列数据统计方法的集合，只不过统计结果用网格形式保存而已。
使用 Data Metrics 方法做网格化一般要进高级选项里面设置相关的参数，最常见的就是统计的方法和搜索的选项。因为我们的目标是做热度图，实则是个点密度图，所以在高级选项里面，统计方法选用 Data Location Statistics 下的 Count 方法。搜索参数简单一些，只设置两个搜索半径为 20 就好，其他采用默认值。
整个网格化参数设置过程，用图片表示如下：
Data Metrics 网格化方法高级选项里面的一些参数设置是环环相扣的，譬如我们选择了 Count 这个统计方法，就必须相应的设置搜索参数。因为 Count 方法的含义是计算搜索区域内的样本（也就是数据点）个数，也就是说，这个方法生成的网格文件里面，每个网格节点值实际上就是以该节点为中心的搜索椭圆内发现的数据个数。如果不设置搜索参数，按 Surfer 的默认设置，是使用全部数据，什么意思呢？意思是， Surfer 默认使用一个巨大的搜索椭圆，大到什么程度？大到这个搜索椭圆以任何一个节点为中心，都能把整个数据集包含在里面。对应到 Count 方法，你就会发现生成的网格没法用啊，因为整个网格里所有节点的值都是同一个数。
也许有人会问，既然是做密度图，为什么不选用 Approximate Density 方法呢？那个看名字不是更适合？实际上是可以。从方法本身而言，Approximate Density 方法等于上述 Count 方法的结果再除以搜索椭圆的面积。两种方法得到的结果，仅仅是在数值上有一个固定的倍数差，在可视化呈现上是完全没有区别的。所以我们选择更简单、更直接、更快的 Count 方法，如此而已。
回头再说说搜索半径。是什么让我们把搜索半径设置为 20 呢？用 5 可以吗？用 30 可以吗？好吧，首先，我们来明确一点，搜索半径里面的数值，所使用的单位与 Z 列一样的。在我们使用的全球地震数据中，在本例当中，我们把 Z 值设为维度列，其单位是度，所以搜索半径也就是经纬度单位度。具体的数值怎么定？通常我们根据自己的研究目的来初步确定，再通过试绘图，找到需要的半径。
本例中，我们通过浏览数据点位图，判断出大致以20°为半径划分搜索区域，会比较具有代表性。再用不同的搜索参数分别网格化，最终确定搜索半径为 20。不信的话，看看用 5 和 30 作为搜索半径会是什么结果：
显然，搜索半径过小，网格节点值过于依赖节点周围的点密度；反过来，搜索半径太大，网格化结果过于宽泛，大片大片数据空白区也被划入热度区，失去代表性。
解决了那么几个乌哩单叨的问题，终于生成了需要的网格文件，赶紧用 Map | New | Contour Map 命令将等值线图画出来。
对于热度图，我们一般不接受默认参数，所以立刻选中等值线图，到属性管理器设置一些参数，主要是 Levels 部分，如下图：
因为我们是用
Count方法做的网格化，所以小于 0 的网格节点值实际上是无效的，因此把Minimum Contour设置为 0 或者 1 就可以了。我们要用大片大片的颜色来表示热区，而 Surfer 默认只显示等值线线条，因此我们把
Fill contours勾选上，跟着在Fill colors挑选一个色谱。自己定做也简单，譬如我就自己用白兰黄红四色定做了一个色谱。有了色块显示就无需等值线及其标注在那里碍眼啦，于是把下面的
Major Contours和Minor Contours的线条都设置为Invisible，把Show labels的勾都去掉。
做完上面的等值线设置，得到下面的图，看起来很粗糙的样子，不满意。
孔夫子说食不厌精烩不厌细，作图也是这个道理，对精细作图要毕生孜孜不倦地追求。对于这种歪模怪样、棱角峥嵘的，在物探上说就是高频成分太多了，要压制。 Surfer 在处理这类问题有光滑和滤波两种选择，既然要消除或压制高频，那我们就用低通滤波好了。
点开 Grid | Filter 命令，直接用默认的 Gaussian (3x3) 低通滤波器及其默认参数，对之前的网格文件做滤波，重新作图，效果是立竿见影：
做完这么光滑圆润的热度图，是不是很有成就感？然并卵！就本例而言，请你告诉我，那个最热（就是地震最高发）地区是在地球上哪个国家哪个地区？不要说慢慢算经纬度然后查表可知。
为了解决这种毫无技术含量然而影响巨大的问题，我们给热度图加一个底图，这样对热度图的判识就简单了。在本例中，我们研究的是全球范围内的地震发生情况，所以加一个简易的世界底图上去。
我用的是从 Golden Software 官网下的 world.gsb 边界文件，但这类全球范围的边界文件网络上随便就能找到，格式也各种各样，够用就好。
既然是边界文件，就用基底图方式，在 Surfer 里面选中刚刚生成的等值线图，然后选择 Map | Add | Base Layer 命令把 world.gsb 边界文件变成底图。因为我的边界文件范围比地震数据范围大，所以会问要不要调整地图限制啊，如下图：
调整个什么！我这是以等值线图为核心的绘图方式，如果按默认调整地图界限，那接下来等值线图就会显得很难看，所以这里要回答“否”。
把底图用灰度色填充（Pattern 设置为 Solid，Foreground color 设置为 30% Black），再把底图放到最后面免得挡住等值线图层（用命令 Arrange | Order Objects | Move To Back），把等值线图层的透明度调到 60%（先选中等值线图层，再到 Property Manager 找 Layer 标签页），这样透过等值线图层就能隐约看到底图了。
为了直观地感受热度图与实测数据之间的关系，可以把原始数据点位图也放上来，用张贴图或者分类张贴图都可以。因为本例的地震数据文件中包含了震级数据，索性做个分类张贴图更加直观和美观。
跟刚才创建基底图一样，添加张贴图层也可能会问你要不要调整地图限制啥的，显然不要，直接点否。然后选中这个分类张贴图，在属性管理器里，把 X、Y、Z 列设置为正确的列即可（本例分别设为经度、纬度和震级列）。再把下面一点的“Show Legend” 勾上就有图例了。
各图层设置完毕，回过头把整张图的其他属性设置一下（譬如比例），把左、下两个坐标的标注设置为经纬度（DMS (lat/long)），在顶坐标设置一个标题。
整张图的完成效果如下：
用脚本绘制热度图
Surfer 具备比较完善的 Automation 功能，自带有使用 VBA 语言的 Scripter 程序，用脚本编程控制 Surfer 自动绘图是相当简单的事情，就是用代码把前面手工完成的事情逐个做了。
而编写脚本代码也很简单，稍微了解一点 VBA 语法，看看 Scripter 的帮助记下几个关键字，然后对着 Surfer 帮助里面的 Automation 部分，把需要的代码拼凑起来就可以。
以下是本例完整代码：
'==================================================================
'CreateHeatMap.bas
'热度图（或密度图）用来可视化特定区域内的密度区、关注区、缓冲区等。
'例如用热度图来显示全球的地震活动区域，虽然全球到处都有地震发生，
'但地震高发区（或易发区）却只是一小部分。
'实例数据可从USGS(http://earthquake.usgs.gov/earthquakes/map/)下载。
'在 Scripter + Golden Software Surfer 13.x 测试通过。
'holz [AT] 21cn.com
'2016-3-22
'==================================================================
Option Explicit '要求变量预先声明
Sub Main
Const COLX = 3 '原始数据中 X 坐标所在列
Const COLY = 2 '原始数据中 Y 坐标所在列
Const COLZ = 2 '原始数据中 Z 坐标所在列
Const COLP = 5 '原始数据中分类标注所在列
'Data Metrics 方法必须设置搜索半径，否则整张图就是一个值。
'如果搜索半径过小，就是各个数据点独立发生作用，反之，如果
'搜索半径设置过大，所得密度图就过于笼统而没有代表性。
'作为示例文件的是全球数据，用经纬度做单位，所以两个半径都设为20。
Const RAD1 = 20
Const RAD2 = 20
Dim SurferApp As Object
Dim PlotDoc As Object
Dim datFile As String '原始数据文件名
Dim grdFile As String '生成的网格文件名
Dim blnFile As String '底图之边界文件名
Dim tmpFile As String '临时文件
Dim bStatus As Boolean '函数执行状态
Dim blkVal As Double '白化值
Dim xLim1 As Double, xLim2 As Double 'X 限制值
Dim yLim1 As Double, yLim2 As Double 'Y 限制值
Dim zLim1 As Double, zLim2 As Double 'Z 限制值
Dim FuncStr As String '数学函数表达式
'已知当前目录下有原始数据文件
datFile = CurDir() & "\2.5_month.csv"
'如果没有，则弹出对话框选择原始数据文件
If Dir(datFile) = "" Then
datFile = GetFilePath(,"csv|dat|txt",,"选择原始数据",4)
End If
'如果仍然未选择任何文件，程序结束
If datFile = "" Then End
Debug.Print "这一切始于： " & datFile
grdFile = Left(datFile,InStrRev(datFile,".")-1)
tmpFile = grdFile & "_tmp.grd"
grdFile = grdFile & ".grd"
'创建一个 Surfer 实例
Set SurferApp = CreateObject("Surfer.Application")
'新增一个空白图形文档
Set PlotDoc=SurferApp.Documents.Add(srfDocPlot)
'==============================================================
'对原始数据进行网格化，网格化方法选择 Data Metrics。
'统计方法选择 Data Location Statistics 下面的 Count。
'启用数据搜索，搜索半径决定每个网格节点涉及的数据个数。
'如需控制结果网格，可设置 xMin, xMax, yMin, yMax,
'NumCols, NumRows, xSize, ySize 等相关参数。
bStatus = SurferApp.GridData(DataFile:=datFile, OutGrid:=grdFile, _
xCol:=COLX, yCol:=COLY, zCol:=COLZ, ShowReport:=False, _
Algorithm:=srfDataMetrics, DataMetric:=srfDMCount, SearchEnable:=True, _
SearchRad1:=RAD1, SearchRad2:=RAD2, SearchAngle:=0, _
SearchNumSectors:=4, SearchMaxData:=64, SearchDataPerSect:=16, _
SearchMinData:=8, SearchMaxEmpty:=3)
Debug.Print "网格化数据" & If(bStatus,"","不") & "成功。"
'对刚生成的网格文件做滤波，使用 Low pass Filter 下的 Gaussian(3x3) 方法。
bStatus = SurferApp.GridFilter(InGrid:=grdFile, OutGrid:=tmpFile, _
Filter:=srfFilterGaussian)
Debug.Print "网格高斯低通滤波" & If(bStatus,"","不") & "成功。"
FuncStr = "C = if(A < 0, 0, A)"
bStatus = SurferApp.gridmath(Function:=FuncStr, InGridA:=tmpFile, OutGridC:=grdFile)
Debug.Print "消除异常网格节点" & If(bStatus,"","不") & "成功。"
'取一些网格参数，用来为后面的绘图服务
Dim Grid As Object
'需要先建立一个 NewGrid 对象，才能塞一个网格进来
Set Grid = SurferApp.NewGrid
With Grid
'只载入网格文件头，免得浪费时间和内存
.LoadFile(FileName:=grdFile, HeaderOnly:=True)
blkVal = .BlankValue
xLim1 = .xMin
xLim2 = .xMax
yLim1 = .yMin
yLim2 = .yMax
zLim1 = .zMin
'对于热度研究，小于 1 就是凉开水吧
If zLim1 < 1 Then zLim1 = 1
zLim2 = .zMax
End With
'一切以等值线图为基准，所以先画等值线
Dim MapFrame As Object
Set MapFrame = PlotDoc.Shapes.AddContourMap(grdFile)
Dim ContourLayer As Object
Set ContourLayer = MapFrame.Overlays(1)
With ContourLayer
'使用简单方法设置等值线
.LevelMethod = SrfConLevelMethodSimple
'设置等值线间隔参数
.SetSimpleLevels(zLim1, zLim2, 2)
'让等值线光滑一些
.SmoothContours = srfConSmoothMed
'将等值线线条设为不可见，因为我们无需等值线
.MajorLine.Style = "Invisible"
.MinorLine.Style = "Invisible"
'也无需等值线标注
.ShowMajorLabels = False
.ShowMinorLabels = False
'要填充等值线必须按一下顺序设置：
'1、启用等值线填充
.FillContours = True
'2、载入一个填充用的色谱
If Dir(CurDir() + "\HeatMaps.clr") <> "" Then
.FillForegroundColorMap.LoadFile(CurDir() + "\HeatMaps.clr")
Else
.FillForegroundColorMap.LoadFile(SurferApp.Path + "\ColorScales\Rainbow.clr")
End If
'3、将填充属性应用到等值线图上
.ApplyFillToLevels(FirstIndex:=1, NumberToSet:=1, NumberToSkip:=0)
'把等值线图层的透明度设置为 60%
.Opacity = 60
End With
blnFile = CurDir() & "\world.gsb"
'如果当前目录有底图的边界文件，就画个底图
If Dir(blnFile) <> "" Then
Dim BaseLayer As Object
'通过添加图层的方式加入底图
Set BaseLayer = PlotDoc.Shapes.AddBaseLayer(MapFrame, blnFile)
'用实色填充闭合边界
BaseLayer.Fill.Pattern = "Solid"
'用 30% 的灰度色进行填充
BaseLayer.Fill.ForeColorRGBA.Color = srfColorBlack30
'底图么，躺到最下面去
BaseLayer.SetZOrder(srfZOToBack)
End If
Dim PostLayer As Object
'用分类张贴图展示实测数据
Set PostLayer = PlotDoc.Shapes.AddClassedPostLayer(MapFrame, datFile, COLX, COLY, COLP)
'对地图容器进行相关设置
Dim Axis As Object
With MapFrame
'用等值线的范围来限定整张图
.SetLimits(xLim1, xLim2, yLim1, yLim2)
'设定XY绘图比例
.xMapPerPU = 10
.yMapPerPU = 10
'对各个坐标轴做相关设置
For Each Axis In .Axes
Select Case Axis.AxisType
'如果是顶坐标，就添加一个标题
Case srfATTop
Axis.Title = "全球地震高发区热度图" & vbCrLf & Date
Axis.TitleFont.Face = "宋体"
Axis.TitleFont.Size = 30
Axis.TitleFont.HAlign = srfTACenter
'坐标刻度不显示
Axis.MajorTickType = srfTickNone
Axis.MinorTickType = srfTickNone
'坐标标注不显示
Axis.ShowLabels = False
Case srfATBottom, srfATLeft
Axis.MajorTickType = srfTickOut
Axis.MinorTickType = srfTickOut
Axis.ShowLabels = True
'可惜目前无法通过 Automation 设置 DMS 形式的标注
'Axis.LabelFormat.Type = srflabdms
Case Else
'坐标刻度不显示
Axis.MajorTickType = srfTickNone
Axis.MinorTickType = srfTickNone
'坐标标注不显示
Axis.ShowLabels = False
End Select
Next
End With
'显示热度图等值线填色色标
ContourLayer.ShowColorScale = True
Dim disColorScale As Object
'获取色标对象
Set disColorScale = ContourLayer.ColorScale
'对色标做相关设置
With disColorScale
If RAD1=RAD2 Then
.Title = "半径 " & RAD1 & " 范围内地震发生次数"
Else
.Title = "长轴 " & IIf(RAD1>RAD2,RAD1,RAD2) & " 短轴 " & _
IIf(RAD1<RAD2,RAD1,RAD2) & " 椭圆范围内地震发生次数"
End If
.TitleFont.Face = "宋体"
.TitleFont.Size = 16
.FirstLabel = 2
.LabelFrequency = 2
.FrameLine.Style = "Invisible"
End With
'显示实测数据分类张贴图图例，因为默认只有方框类型的图例，
'想要实现水平横向分布的图例展示形式，只好自己解决。
Dim x As Double, y As Double, w As Double
Dim LegendTitle As Object
Dim xInc1 As Double, xInc2 As Double
xInc1 = 0.5 '图例符号和图例说明之间的间距为 0.5cm
xInc2 = 1.0 '图例之间的间距为 1.0cm
x = MapFrame.Left + MapFrame.Width * 0.5 '地图中心的 X 坐标
y = MapFrame.Top - MapFrame.Height - 1.0 '地图下方 1cm 处的 Y 坐标
'先写图例的标题
Set LegendTitle = PlotDoc.Shapes.AddText(x,y,"实测数据分类图例：")
With LegendTitle
.Font.Face = "宋体" '图例标题用中文字体
.Font.Size = 16 '文本的大小
.Font.HAlign = srfTALeft '水平靠左对齐
.Font.VAlign = srfTAVCenter '垂直居中对齐
w = .Width + xInc1 '统计整个图例内所有对象的宽度总和
End With
Dim NumBins As Integer
Dim LegendSym() As Object, LegendLab() As Object '动态数组预定义
Dim i As Integer
Dim txtStr As String
NumBins = PostLayer.NumClasses '分类张贴图的符号组数
ReDim LegendSym(NumBins), LegendLab(NumBins) '给动态数组规定大小
'把每个组的相关参数取出来，做成图例
For i = 1 To NumBins
'在图上画一个图例符号
Set LegendSym(i) = PlotDoc.Shapes.AddSymbol(x,y)
With LegendSym(i).Marker
'他的符号集就是对应分组使用的符号集
.Set = PostLayer.BinSymbol(i).Set
'他所用的符号就是对应分组使用的符号
.Index = PostLayer.BinSymbol(i).Index
'他的大小就是对应分组的符号大小
.Size = PostLayer.BinSymbol(i).Size
'他的颜色就是对应分组的符号颜色
.LineColor = PostLayer.BinSymbol(i).LineColor
.FillColor = PostLayer.BinSymbol(i).FillColor
End With
w = w + LegendSym(i).Width + xInc1
'图例标注用“分组下限 - 分组上限”的形式，可随意修改
txtStr = PostLayer.BinLowerLimit(i) & " - " & PostLayer.BinUpperLimit(i)
'把图例标注画到图上
Set LegendLab(i) = PlotDoc.Shapes.AddText(x,y,txtStr)
'设置图例标注的字体、大小、对齐等属性
With LegendLab(i)
.Font.Size = 12
.Font.HAlign = srfTALeft
.Font.VAlign = srfTAVCenter
w = w + .Width + xInc2
End With
Next
'图例生成之后，就要按需要摆放，本例是水平横向摆放。
'由于其位置参照的是地图本身，所以应该在设置好整张图的
'比例大小之后再摆放图例，顺序错了，图例位置就不对。
'下列代码实现的是在图形下方水平居中横向布设图例。
Dim x0 As Double, w0 As Double
x0 = MapFrame.Left '取出地图容器的左坐标
w0 = MapFrame.Width '取出地图容器的宽度
x = x0 + (w0 - w) * xInc1
LegendTitle.Left = x
x = x + LegendTitle.Width + xInc1
For i = 1 To NumBins
LegendSym(i).Left = x
x = x + LegendSym(i).Width + xInc1
LegendLab(i).Left = x
x = x + LegendLab(i).Width + xInc2
Next
'如果用软件自带的图例功能，只需要一下几行代码即可
'告诉张贴图层要显示图例
'PostLayer.ShowLegend = True
'Dim pstLegend As Object
'获取图例对象
'Set pstLegend = PostLayer.Legend
'对图例做各种设置
'With pstLegend
' .Title = "实测数据图例"
' .TitleFont.Face = "宋体"
' .TitleFont.Size = 16
'End With
'让图形缩放到适应窗口
SurferApp.ActiveWindow.Zoom(srfZoomFitToWindow)
'让 Surfer 主窗口可见
SurferApp.Visible=True
Debug.Print "以下两个过程文件可删除：" & vbCrLf & _
" --> " & tmpFile & vbCrLf & _
" --> " & grdFile
End Sub
用这个脚本绘制的热度图如下：
小结
手工绘图和自动化绘图各有优缺点。手工绘图在于每一步几乎都是所见即所得，自动化绘图在于只要写好代码随时一键成图。
绘图步骤很多或者要绘制大量同类图件时，编写脚本是提高效率避免出错的不二之选。
自动化绘图也不是万能，必要时还是要做后期的手工调整，例如用脚本就无法设置坐标标注为度分秒形式。