Средние значения набора точечных данных для набора данных сетки

Я относительно новичок в ggplot, поэтому, пожалуйста, простите меня, если некоторые из моих проблем действительно просты или вообще не решаемы.

Я пытаюсь создать «тепловую карту» страны, в которой фигура заполняется непрерывно. Кроме того, у меня есть форма страны как .RData. Чтобы преобразовать мою Данные 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 Netherlands


person Marco    schedule 19.12.2011    source источник
comment
Приведите воспроизводимый пример.   -  person Paul Hiemstra    schedule 19.12.2011
comment
У меня такое чувство, что вы ищете алгоритм интерполяции, см. Мой пост ниже, чтобы увидеть пример использования кригинга (геостатистики).   -  person Paul Hiemstra    schedule 19.12.2011
comment
Отлично, что вы опубликовали решение, +1. Единственное, на что я хотел бы обратить внимание, это отсутствие библиотеки (масштабов) для отключенной функции.   -  person Eduardo    schedule 02.06.2013


Ответы (2)


Я думаю, что вы хотите что-то в этом роде. Я предсказываю, что это домашнее пиво будет ужасно неэффективным для больших наборов данных, но оно работает с небольшим примером набора данных. Я бы посмотрел на плотность ядра и, возможно, на пакет 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()

что приводит к следующему изображению:

введите описание изображения здесь

Альтернативный подход - использовать алгоритм интерполяции. Один из примеров - кригинг. Это довольно просто, используя пакет 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()

что приводит к следующему изображению:

введите описание изображения здесь

Однако, не зная, какова ваша цель с этой картой, мне трудно понять, что именно вы хотите.

person Paul Hiemstra    schedule 19.12.2011

Это слайд-шоу предлагает другой подход - см. стр. 18 для описания подхода и стр. 21, чтобы узнать, как выглядели результаты для создателя слайдов.

Однако обратите внимание, что средство создания слайдов использовало пакет sp и функцию spplot, а не ggplot2 и его функции построения графиков.

person prabhasp    schedule 20.12.2011
comment
... кроме того, spplot использует решетку под капотом. - person Paul Hiemstra; 20.12.2011