在地图上绘制插值数据

Question

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

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

首先，谁能评论我是否正确实现了autoKrige() 函数？

其次，我在绘制数据和覆盖该地区的地图时遇到了麻烦。或者，我可以指定插值网格来反映海湾的边界（如建议的here）吗？关于我如何做到这一点以及我可以从哪里获得这些信息的任何想法？将网格提供给 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")

Answer 1

我对你的帖子有一些评论：

使用克里金法

我看到您正在使用地统计学来构建热图。您还可以考虑其他插值技术，例如样条曲线（例如，字段包中的薄板样条曲线）。这些对数据的假设更少（例如平稳性），并且还可以很好地可视化您的数据。如果您将其发送到期刊，那么减少假设数量可能会有所帮助，这样您就无需向审稿人解释了。如果需要，您还可以比较一些插值技术，请参阅a report I wrote 了解一些提示。

数据投影

我看到您正在使用经纬度坐标进行克里金法。 Edzer Pebesma（gstat 的作者）remarked that 没有适合经纬度坐标的变异函数模型。这是因为在 lat lon 中，距离不是直线（即Euclidean），而是在一个球体上（即Great circle distances）。没有对球坐标有效的协方差函数（或变异函数模型）。我建议在使用自动映射之前使用 rgdal 包中的 spTransform 投影它们。

rgdal 包使用proj4 projection library 来执行计算。要投影您的数据，您首先需要定义它的投影：

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

上面表达式右侧的 proj4 字符串定义了投影的类型 (+proj)、使用的省略号 (+ellps) 和基准面 (+datum)。要理解这些术语的含义，您必须将地球想象成一个土豆。地球不是完美的球形，这是由椭圆定义的。地球也不是一个完美的椭球体，但表面更不规则。这种不规则性由基准定义。另见this article on Wikipedia。

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

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

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

用 ggplot2 绘图

要在 ggplot 中添加多边形或折线，请查看coord_map 的文档。这包括使用maps 包绘制国家边界的示例。如果您需要为您的研究区域加载例如 shapefile，您可以使用rgdal 来完成。请记住 ggplot2 适用于 data.frame，而不是 SpatialPolygons。您可以使用以下方法将SpatialPolygons 转换为data.frame：

poly_df = fortify(poly_Spatial)

另见this function 我创建用于绘制空间网格。它直接在 SpatialGrids/Pixels 上工作。请注意，您需要从该存储库 (continuousToDiscrete) 中获取一两个附加文件。

创建插值网格

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

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

如果您将插值地图用于后续分析，则获得正确的分辨率很重要。如果您将地图纯粹用于可视化目的，那么这一点就不那么重要了。但请注意，在这两种情况下，过高的分辨率都可能会误导您的预测准确性，而过低的分辨率不会对数据产生公正的影响。