Valores promedio de un conjunto de datos de punto a un conjunto de datos de cuadrícula

Question

Soy relativamente nuevo en ggplot, así que, por favor, discúlpeme si algunos de mis problemas son realmente simples o no se pueden resolver.

Lo que intento hacer es generar un "Mapa de calor" de un país donde el relleno de la forma es continuo. Además tengo la forma del país como .RData. Usé hadley wickham's script para transformar mis datos de SpatialPolygon en un marco de datos. Los datos largos y lat de mi marco de datos ahora se ven así

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 

Estos datos long/lat dibujan el contorno de Alemania. El resto del marco de datos se omite aquí porque creo que no es necesario. También tengo un segundo marco de datos de valores para ciertos puntos largos/Lat. Esto se parece a esto

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 

Lo que me gustaría hacer ahora, es el color de cada punto del mapa de acuerdo con un valor promedio de todos los puntos fijos que se encuentran dentro de una cierta distancia de ese punto. De esa forma obtendría una coloración (casi) continua de todo el mapa del país. Lo que tengo hasta ahora es el mapa del país trazó con ggplot2

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

Mi primera idea era generar puntos que están dentro del mapa que se ha dibujado y luego calcular el valor de cada punto my_generated_point generada como tal

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) 

No he encontrado una manera de generar estos puntos. Y como con todo el problema, ni siquiera sé si es posible resolverlo de esta manera. Mi pregunta ahora es si existe una forma de generar estos puntos y/o si hay otra manera de llegar a una solución.

Solución

Gracias a Paul casi me tengo lo que quería. Aquí hay un ejemplo con datos de muestra para los Países Bajos.

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

Creo que lo que quiere es algo a lo largo de estas líneas. Predigo que este homebrew va a ser terriblemente ineficiente para grandes conjuntos de datos, pero funciona en un pequeño conjunto de datos de ejemplo. Me gustaría ver las densidades del núcleo y tal vez el paquete raster. Pero quizás esto te convenga bien ...

El siguiente fragmento de código calcula el valor medio de la concentración de cadmio de una grilla de puntos que se superpone al conjunto de datos de puntos original. Solo se consideran puntos a menos de 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() 

que conduce a la siguiente imagen:

Un enfoque alternativo es utilizar un algoritmo de interpolación. Un ejemplo es kriging.Esto es muy fácil usando el paquete automapa (:) detectar la auto-promoción, escribí el paquete):

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

que conduce a la siguiente imagen:

Sin embargo, sin el conocimiento como para cuál es tu objetivo con este mapa, es difícil para mí ver exactamente lo que quieres.

Answer 2

Este slideshow ofrece otro enfoque: consulte page 18 para obtener una descripción de la aproximación y page 21 para ver los resultados del deslizador.

Sin embargo, tenga en cuenta que el fabricante de diapositivas usó el paquete sp y la función spplot en lugar de ggplot2 y sus funciones de trazado.