2011-12-19 110 views
8

我對ggplot比較新,所以請原諒我,如果我的一些問題很簡單或根本無法解決。點數據集到網格數據集的平均值

我想要做的是生成一個國家的「熱圖」,形狀的填充是連續的。此外,我的國家形狀爲.RData。我使用hadley wickham's script將SpatialPolygon數據轉換爲數據框。我的數據幀的長和經度數據現在看起來像這樣

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 

這個long/lat數據繪製了德國的輪廓。數據框的其餘部分在這裏省略,因爲我認爲它不是必需的。對於特定的長/經點,我還有第二個數據框。這看起來像這樣

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 Netherlands

+0

請做一個可重現的例子。 – 2011-12-19 15:17:08

+0

我有一種感覺,你正在尋找一種插值算法,請參閱下面我的帖子,以克里格(地質統計學)爲例。 – 2011-12-19 15:53:34

+0

很棒,你已經發布瞭解決方案,+1。我唯一要指出的是它缺少「靜音」功能的庫(縮放)。 – Eduardo 2013-06-02 19:36:16

回答

15

我想你想要的是沿着這些路線的東西。我預測這種自制軟件對於大型數據集的效率會非常低,但它可以用於一個小的示例數據集。我會研究內核密度,也許是raster包。但也許這適合你...

下面的代碼片段計算覆蓋原始點數據集點的網格鎘濃度的平均值。只考慮接近1000米的點。

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() 

這導致以下圖片:

enter image description here

的另一種方法是使用內插算法。克里格就是一個例子。這是使用自動地圖包很容易(臨場自我宣傳:),我寫的包):

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() 

這導致了下面的圖片:

enter image description here

然而,如果沒有知識到這張地圖你的目標是什麼,我很難看到你想要的東西。

2

slideshow提供了另一種方法 - 有關該方法的描述的page 18和有關幻燈片製造商的結果看起來如何的page 21

但請注意,幻燈片製作工具使用sp包和spplot函數,而不是ggplot2及其繪圖函數。

+0

...此外,spplot使用引擎蓋下的格子。 – 2011-12-20 10:35:36