Giá trị trung bình của tập dữ liệu điểm cho một tập dữ liệu lưới

Question

Tôi tương đối mới đối với ggplot, vì vậy hãy tha thứ cho tôi nếu một số vấn đề của tôi thực sự đơn giản hoặc không thể giải quyết được.

Điều tôi đang cố gắng tạo ra là "Bản đồ nhiệt" của một quốc gia nơi việc điền hình dạng liên tục. Hơn nữa tôi có hình dạng của đất nước là .RData. Tôi đã sử dụng hadley wickham's script để chuyển đổi dữ liệu SpatialPolygon của tôi thành một khung dữ liệu. Dữ liệu dài và lat của khung dữ liệu của tôi bây giờ trông giống như thế này

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 

Dữ liệu dài/lat này vẽ phác thảo của Đức. Phần còn lại của khung dữ liệu bị bỏ qua ở đây vì tôi nghĩ rằng nó không cần thiết. Tôi cũng có một khung dữ liệu thứ hai của các giá trị cho một số điểm dài/lat nhất định. Điều này trông giống như thế này

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 

Điều tôi muốn làm bây giờ là tô màu mỗi điểm trên bản đồ theo một khoảng cách nhất định. Bằng cách đó tôi sẽ nhận được một màu (liên tục) liên tục của toàn bộ bản đồ của đất nước. Những gì tôi có cho đến nay là bản đồ của nước này âm mưu với ggplot2

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

Idea đầu tiên của tôi là để tạo ra điểm nằm trong bản đồ đã được vẽ và sau đó tính toán giá trị cho tất cả các điểm tạo my_generated_point như vậy

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) 

Tôi chưa tìm được cách tạo ra các điểm này. Và như với toàn bộ vấn đề, tôi thậm chí không biết nếu nó có thể giải quyết theo cách này. Câu hỏi của tôi bây giờ là nếu có một cách để tạo ra những điểm này và/hoặc nếu có một cách khác để đi đến một giải pháp.

Giải pháp

Nhờ Paul Tôi gần như nhận được những gì tôi muốn. Đây là một ví dụ với dữ liệu mẫu cho Hà Lan.

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

Answer 1

Tôi nghĩ rằng những gì bạn muốn là một cái gì đó dọc theo những dòng này. Tôi dự đoán rằng homebrew này sẽ vô cùng hiệu quả đối với các tập dữ liệu lớn, nhưng nó hoạt động trên một tập dữ liệu mẫu nhỏ. Tôi sẽ xem xét mật độ hạt nhân và có thể là gói raster. Nhưng có thể điều này phù hợp với bạn ...

Đoạn mã sau đây tính toán giá trị trung bình của nồng độ cadmium của một mạng lưới các điểm phủ lên tập dữ liệu điểm gốc. Chỉ những điểm gần hơn 1000 m mới được xem xét.

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

dẫn đến các hình ảnh sau đây:

Cách tiếp cận khác là sử dụng một thuật toán nội suy. Một ví dụ là kriging.Điều này khá dễ dàng bằng cách sử dụng gói AutoMap (tại chỗ :) tự quảng bá, tôi đã viết các gói):

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

dẫn đến hình ảnh sau:

Tuy nhiên, không có kiến ​​thức như với mục tiêu của bạn là gì với bản đồ này, thật khó cho tôi để xem những gì bạn muốn chính xác.

Answer 2

slideshow cung cấp cách tiếp cận khác - xem page 18 để biết mô tả về cách tiếp cận và page 21 để biết kết quả trông như thế nào đối với trình tạo bản chiếu.

Lưu ý rằng trình tạo bản chiếu đã sử dụng gói sp và hàm spplot thay vì ggplot2 và các chức năng vẽ bản đồ của nó.