jsl*_*che 24 maps r ggplot2 automap spatial-interpolation
我有在美国切萨皮克湾不同地点拍摄的物种丰富度的调查数据,我想以图形方式将数据显示为"热图".
我有一个lat/long坐标和丰富度值的数据框,我将其转换为a SpatialPointsDataFrame并使用autoKrige()automap包中的函数生成插值.
首先,任何人都可以评论我是否正确实现了该autoKrige()功能？
其次,我无法绘制数据并覆盖该地区的地图.或者,我可以指定插值网格来反映Bay的边界(如此处所示)吗？关于我如何做到这一点以及我可能从哪里得到这些信息的任何想法？提供网格autoKrige()看起来很容易.
编辑:感谢Paul的超级有用的帖子!这就是我现在拥有的.无法让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")
Pau*_*tra 47
我对你的帖子有很多评论:
我看到您正在使用地统计学来构建热图.您还可以考虑其他插值技术,例如样条线(例如,字段包中的薄板样条).这些对数据的假设较少(例如平稳性),也可以很好地可视化您的数据.如果您将其发送到期刊,假设数量的减少可能会有所帮助,那么您就可以向审稿人解释.如果需要,您还可以比较一些插值技术,请参阅我为一些提示撰写的报告.
我看到你正在使用lat long坐标进行克里金法.Edzer Pebesma(作者gstat)评论说,没有适合拉特伦坐标的变差函数模型.这是因为在纬度方向上距离不是直的(即欧几里得),而是在一个球体上(即大圆距离).没有协方差函数(或变差函数模型)对球坐标有效.我建议在使用automap之前spTransform从rgdal包中使用它们进行投影.
rgdal包使用proj4投影库来执行计算.要投影数据,首先需要定义其投影:
proj4string(df) = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
上面表达式右侧的proj4字符串定义了投影(+proj)的类型,使用的椭圆(+ellps)和基准(+datum).要理解这些术语的含义,你必须将地球想象成一个土豆.地球不是完美的球形,这是由椭圆形定义的.地球也不是完美的椭圆体,但表面更不规则.这种不规则性由基准定义.另见维基百科上的这篇文章.
定义投影后,您可以使用spTransform:
project_df = spTransform(df, CRS("+proj= etcetc"))
其中CRS("+ proj etc")定义了目标投影.哪种投影取决于您的地理位置和学习区域的大小.
要将多边形或折线添加到ggplot,请查看文档coord_map.这包括使用maps包绘制国家边界的示例.如果您需要为您的研究区域加载例如shapefile,您可以使用rgdal.请记住,它ggplot2适用于data.frame,而不是SpatialPolygons.您可以转换SpatialPolygons为data.frame使用:
poly_df = fortify(poly_Spatial)
另请参阅我创建的用于绘制空间网格的函数.它直接适用于SpatialGrids/Pixels.请注意,您需要从该存储库中获取一个或两个附加文件(continuousToDiscrete).
当没有指定时,我创建了automap来生成输出网格.这是通过在数据点周围创建一个凸包,并在其中采样5000点来完成的.预测区域的边界以及在其中采样的点的数量(以及因此的分辨率)是非常随意的.对于特定应用,预测区域的形状可以从多边形导出,spsample用于采样多边形内的点.取样的数量和分辨率取决于两件事:
如果您使用插值地图进行后续分析,那么获得正确的分辨率非常重要.如果您纯粹将地图用于视觉目的,那么这一点就不那么重要了.但请注意,在这两种情况下,过高的分辨率可能会误导预测的准确性,而分辨率太低则无法确定数据的正确性.