24

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

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

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

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

1 回答 1

47

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

使用克里金法

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

数据投影

我看到您正在使用经纬度坐标进行克里金法。Edzer Pebesma（作者gstat）评论说，没有适合经纬度坐标的变异函数模型。这是因为在纬度中，距离不是直线（即欧几里得），而是在一个球体上（即大圆距离）。没有对球坐标有效的协方差函数（或变异函数模型）。我建议在使用 automap 之前spTransform从包中投影它们。rgdal

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

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

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

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

创建插值网格

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

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

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

于 2012-04-07T22:18:41.723 回答