8

我对ggplot比较陌生，所以如果我的一些问题真的很简单或根本无法解决，请原谅我。

我想做的是生成一个国家的“热图”，其中形状的填充是连续的。此外，我有国家的形状.RData。我使用hadley wickham 的脚本将我的 SpatialPolygon 数据转换为数据框。我的数据框的 long 和 lat 数据现在看起来像这样

head(my_df)
long        lat         group
6.527187    51.87055    0.1 
6.531768    51.87206    0.1
6.541202    51.87656    0.1
6.553331    51.88271    0.1

这个长/纬度数据勾勒出德国的轮廓。数据框的其余部分在此省略，因为我认为不需要。对于某些长/纬度点，我还有第二个值数据框。这看起来像这样

my_fixed_points
long        lat         value
12.817      48.917      0.04 
8.533       52.017      0.034
8.683       50.117      0.02
7.217       49.483      0.0542

我现在想做的是根据位于该点一定距离内的所有固定点的平均值为地图的每个点着色。这样我就可以（几乎）连续地为整个国家的地图着色。到目前为止，我所拥有的是用 ggplot2 绘制的国家地图

ggplot(my_df,aes(long,lat)) + geom_polygon(aes(group=group), fill="white") + 
geom_path(color="white",aes(group=group)) + coord_equal()

我的第一个想法是生成位于已绘制地图内的点，然后my_generated_point像这样计算每个生成点的值

value_vector <- subset(my_fixed_points, 
  spDistsN1(cbind(my_fixed_points$long, my_fixed_points$lat),  
  c(my_generated_point$long, my_generated_point$lat), longlat=TRUE) < 50, 
  select = value)
point_value <- mean(value_vector)

我还没有找到产生这些点的方法。和整个问题一样，我什至不知道是否可以通过这种方式解决。我现在的问题是是否存在产生这些点的方法和/或是否有另一种方法来解决问题。

解决方案

多亏了保罗，我几乎得到了我想要的。以下是荷兰的示例数据示例。

library(ggplot2)
library(sp)
library(automap)
library(rgdal)
library(scales)

#get the spatial data for the Netherlands
con <- url("http://gadm.org/data/rda/NLD_adm0.RData")
print(load(con))
close(con)

#transform them into the right format for autoKrige
gadm_t <- spTransform(gadm, CRS=CRS("+proj=merc +ellps=WGS84"))

#generate some random values that serve as fixed points
value_points <- spsample(gadm_t, type="stratified", n = 200)
values <- data.frame(value = rnorm(dim(coordinates(value_points))[1], 0 ,1))
value_df <- SpatialPointsDataFrame(value_points, values)

#generate a grid that can be estimated from the fixed points
grd = spsample(gadm_t, type = "regular", n = 4000)
kr <- autoKrige(value~1, value_df, grd)
dat = as.data.frame(kr$krige_output)

#draw the generated grid with the underlying map
ggplot(gadm_t,aes(long,lat)) + geom_polygon(aes(group=group), fill="white") + geom_path(color="white",aes(group=group)) + coord_equal() + 
geom_tile(aes(x = x1, y = x2, fill = var1.pred), data = dat) + scale_fill_continuous(low = "white", high = muted("orange"), name = "value")

autoKrige 荷兰

4

2 回答 2

15

我认为你想要的是这些方面的东西。我预测这个自制程序对于大型数据集将非常低效，但它适用于小型示例数据集。我会研究内核密度，也许还有raster包。但也许这很适合你...

以下代码片段计算覆盖原始点数据集的点网格的镉浓度平均值。仅考虑接近 1000 m 的点。

library(sp)
library(ggplot2)
loadMeuse()

# Generate a grid to sample on
bb = bbox(meuse)
grd = spsample(meuse, type = "regular", n = 4000)
# Come up with mean cadmium value
# of all points < 1000m.
mn_value = sapply(1:length(grd), function(pt) {
  d = spDistsN1(meuse, grd[pt,])
  return(mean(meuse[d < 1000,]$cadmium))
})

# Make a new object
dat = data.frame(coordinates(grd), mn_value)
ggplot(aes(x = x1, y = x2, fill = mn_value), data = dat) + 
   geom_tile() + 
   scale_fill_continuous(low = "white", high = muted("blue")) + 
   coord_equal()

这导致以下图像：

在此处输入图像描述

另一种方法是使用插值算法。一个例子是克里金法。这很容易使用 automap 包（发现自我提升 :)，我编写了包）：

library(automap)
kr = autoKrige(cadmium~1, meuse, meuse.grid)
dat = as.data.frame(kr$krige_output)

ggplot(aes(x = x, y = y, fill = var1.pred), data = dat) + 
   geom_tile() + 
   scale_fill_continuous(low = "white", high = muted("blue")) + 
   coord_equal()

这导致以下图像：

在此处输入图像描述

但是，如果不知道您对这张地图的目标是什么，我很难确切地看到您想要什么。

于 2011-12-19T15:41:03.703 回答
2

该幻灯片提供了另一种方法——参见第 18 页了解该方法的说明，参见第21 页了解幻灯片制作者的结果。

但是请注意，幻灯片制作器使用了 sp 包和spplot函数，而不是 ggplot2 及其绘图函数。

于 2011-12-20T07:43:33.877 回答