我有在美国切萨皮克湾不同地点采集的物种丰富度调查数据，我想以图形方式将这些数据呈现为“热图”。

我有一个纬度/经度坐标和丰富度值的数据框，我将其转换为SpatialPointsDataFrame并使用了autoKrige()automap 包中的函数用于生成插值。

首先，任何人都可以评论我是否正确实施了autoKrige()功能？

其次，我在绘制数据和覆盖该地区的地图时遇到困难。或者，我可以指定插值网格来反映海湾的边界（如建议的那样）here https://stackoverflow.com/questions/8563334/average-values-of-a-point-dataset-to-a-grid-dataset）？关于如何做到这一点以及从哪里可以获得这些信息有什么想法吗？向电网供电autoKrige()看起来很容易。

编辑：感谢保罗的超级有用的帖子！这是我现在所拥有的。让 ggplot 接受插值数据和地图投影时遇到问题：

require(rgdal)
require(automap)
#Generate lat/long coordinates and richness data
set.seed(6)
df=data.frame(
  lat=sample(seq(36.9,39.3,by=0.01),100,rep=T),
  long=sample(seq(-76.5,-76,by=0.01),100,rep=T),
  fd=runif(10,0,10))
initial.df=df

#Convert dataframe into SpatialPointsDataFrame
coordinates(df)=~long+lat

#Project latlong coordinates onto an ellipse
proj4string(df)="+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
#+proj = the type of projection (lat/long)
#+ellps and +datum = the irregularity in the ellipse represented by planet earth

#Transform the projection into Euclidean distances
project_df=spTransform(df, CRS("+proj=merc +zone=18s +ellps=WGS84 +datum=WGS84")) #projInfo(type="proj")

#Perform the interpolation using kriging
kr=autoKrige(fd~1,project_df)
#Extract the output and convert to dataframe for easy plotting with ggplot2
kr.output=as.data.frame(kr$krige_output)
#Plot the output
#Load the map data for the Chesapeake Bay
cb=data.frame(map("state",xlim=range(initial.df$long),ylim=range(initial.df$lat),plot=F)[c("x","y")])

ggplot()+
  geom_tile(data=kr.output,aes(x=x1,y=x2,fill=var1.pred))+  
  geom_path(data=cb,aes(x=x,y=y))+
  coord_map(projection="mercator")

我对你的帖子有几点评论：

使用克里格法

我看到您正在使用地统计学来构建热图。您还可以考虑其他插值技术，例如样条线（例如 fields 包中的薄板样条线）。这些对数据的假设（例如平稳性）较少，并且还可以很好地可视化您的数据。如果您将其发送给期刊，则假设数量的减少可能会有所帮助，这样您就无需向审稿人解释。如果需要，您还可以比较一些插值技术，请参阅我写的一份报告 http://www.numbertheory.nl/files/report_evap.pdf一些提示。

数据投影

我发现您正在使用经纬度坐标进行克里格法。埃德泽·佩贝斯玛（Edzer Pebesma）（作者gstat) 指出 http://r-sig-geo.2731867.n2.nabble.com/Great-Circle-distances-in-Automap-Gstat-td6863940.html没有适合经纬度坐标的变异函数模型。这是因为在经纬度中，距离不是直的（即欧几里得 http://en.wikipedia.org/wiki/Euclidean_distance），但是在一个球体上，（即大圆距离 http://en.wikipedia.org/wiki/Great-circle_distance）。没有对球坐标有效的协方差函数（或变差函数模型）。我建议使用投影它们spTransform来自rgdal使用 automap 之前先打包。

rgdal 包使用proj4 投影库 http://trac.osgeo.org/proj/执行计算。要投影数据，您首先需要定义其投影：

proj4string(df) = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

上面表达式右侧的 proj4 字符串定义了投影的类型 (+proj），使用的椭圆（+ellps）和数据（+datum）。要理解这些术语的含义，您必须将地球想象成一个土豆。地球不是完美的球形，这是由椭圆定义的。地球也不是一个完美的椭球体，但其表面更加不规则。这种不规则性是由数据定义的。也可以看看维基百科上的这篇文章 http://en.wikipedia.org/wiki/Map_projection.

定义投影后，您可以使用spTransform:

project_df = spTransform(df, CRS("+proj= etcetc"))

其中 CRS("+proj etc") 定义目标投影。哪种投影合适取决于您的地理位置和研究区域的大小。

使用 ggplot2 绘图

要向 ggplot 添加多边形或多段线，请查看以下文档coord_map。这包括使用的示例maps绘制国家边界的包。如果您需要为您的研究区域加载例如形状文件，您可以使用rgdal。请记住这一点ggplot2适用于 data.frame，不适用于SpatialPolygons。你可以转变SpatialPolygons to data.frame using:

poly_df = fortify(poly_Spatial)

也可以看看这个功能 https://bitbucket.org/paulhiemstra/paultools/src/ac65b18511f0/plotSpatialGrids.r我创建是为了绘制空间网格。它直接在 SpatialGrids/Pixels 上工作。请注意，您需要从该存储库中获取一两个附加文件（连续到离散 https://bitbucket.org/paulhiemstra/paultools/src/ac65b18511f0/continuousToDiscrete.r).

创建插值网格

当没有指定时，我创建了 automap 来生成输出网格。这是通过在数据点周围创建凸包并在其中采样 5000 个点来完成的。预测区域的边界以及其中采样的点数（以及分辨率）是相当任意的。对于特定应用，预测区域的形状可以从多边形导出，使用spsample对多边形内的点进行采样。采样多少个点以及分辨率取决于两件事：

• 例如，如果您的数据非常平滑，那么与这种平滑度相比，将分辨率提高得非常高就没有多大意义。或者，如果您的数据具有许多小规模结构，则需要高分辨率。当然，只有当您有支持这种高分辨率的观察结果时，这才有可能。
• 数据的密度。如果您的数据更密集，您可以提高分辨率。

如果您使用插值图进行后续分析，获得正确的分辨率非常重要。如果您纯粹将地​​图用于可视化目的，那么这一点就不那么重要了。但请注意，在这两种情况下，太高的分辨率可能会误导预测的准确性，而太低的分辨率则无法公正地对待数据。