Проблемы с использованием stat_countour () на карте (ggplot2. R)

Question

У меня есть фрейм данных с местоположениями (широта + долгота) в Евразии, и я хотел бы использовать третью переменную для создания контурной карты. Фрейм данных (datgeo) выглядит следующим образом:

lat       long      date_BP
56.28000  25.13000  7429.833
40.31992  29.45311  8048.077
50.41027  14.07460  4200.000
50.12175  14.45695  4484.600
58.74000  -2.91600  4913.444
44.53000  22.05000  8200.333
50.09000  74.44000  3707.125
34.75146  72.40194  2834.625
...

Я использую ggplot2 для создания карты Евразии. И я попытался использовать date_BP в качестве оси z для создания карты контуров

library (ggplot2)
library(rnaturalearth)
library(rnaturalearthdata)
library(ggspatial)
world <- ne_countries(scale = "medium", returnclass = "sf")
ggplot(data = world) +
  geom_sf() +
  geom_point(data = datgeo, aes(x = long, y = lat)) +
  stat_contour(geom="polygon",data=datgeo, aes(x=long,y=lat,z=date_BP,fill=..level..)) +
  scale_fill_gradient(low="blue",high="red")+
  coord_sf(xlim = c(-12.3, 110), ylim = c(70, 30), expand = FALSE)

Однако скрипт не работал вообще, и я не мог получить карту контуров. Это карта, которую генерирует вышеприведенный скрипт:

Почему этот скрипт не представляет никаких контурных линий (изобар)?

Answer 1

Я думаю, что ваша проблема в том, что stat_contour не работает, потому что ему нужна полная сетка. Я нашел статью этого блога, в которой объясняется, как бороться с этой проблемой: https://www.r-statistics.com/2016/07/using-2d-contour-plots-within-ggplot2-to-visualize-relationships-between-three-variables/

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

Сначала необходимо создать прогнозируемую модель на основе вашего ограниченного набора данных "datgeo".

data_geo_loess <- loess(date_BP ~lat+long, data = datgeo)

Затем вы можете создать сетку значений с широтными и длинными значениями:

lat_grid <- seq(min(datgeo$lat),max(datgeo$lat),0.1)
long_grid <- seq(min(datgeo$long), max(datgeo$long),0.1)
data_grid <- expand.grid(lat = lat_grid, long = long_grid)

Теперь вы можете использовать модель loess для расчета теоретических значений на основе date_BP. для всех значений lat и long, которые вы сгенерировали, и мы изменим форму, чтобы получить подходящий фрейм данных для ggplot2:

geo_fit <- predict(data_geo_loess, newdata = data_grid)

library(reshape2)
geo_fit <- melt(geo_fit, id.vars = c("lat","long"), measure.vars = "date_BP")

library(stringr)
geo_fit$lat <- as.numeric(str_sub(geo_fit$lat, str_locate(geo_fit$lat, "=")[1,1] + 1))
geo_fit$long <- as.numeric(str_sub(geo_fit$long, str_locate(geo_fit$long, "=")[1,1] + 1))

> head(geo_fit)
       lat   long    value
1 34.75146 -2.916 24170.02
2 34.85146 -2.916 24290.79
3 34.95146 -2.916 24381.19
4 35.05146 -2.916 24442.12
5 35.15146 -2.916 24474.53
6 35.25146 -2.916 24479.34

Наконец, вы можете получить свой график, выполнив:

library(sf)
library(sp)
library(maps)
library(rnaturalearth)

ggplot(data = world) +
  geom_sf() +
  coord_sf(xlim = c(-12.3, 110), ylim = c(70, 30), expand = FALSE) +
  stat_contour(geom="polygon",
               inherit.aes = FALSE, 
               data=geo_fit, alpha = 0.5, fill = NA,
               aes(x=long,y=lat,z=value, color=..level..)) +
  geom_point(data = datgeo, aes(x = long, y = lat)) +
  scale_color_gradient(low="blue",high="red")

Выглядит ли то, что вы ожидаете?

---

Примечание: loess модель выдаст некоторые предупреждения (по крайней мере, в мой случай) потому что слишком мало наблюдений, чтобы построить надежную модель. Таким образом, вам придется посмотреть с вашими реальными и более полными данными, если они работают.

Примечание: альтернативным решением будет использование stat_density_2d, но вы не можете использовать трехмерное значение.

---

Воспроизводимый пример

structure(list(lat = c(56.28, 40.31992, 50.41027, 50.12175, 58.74, 
44.53, 50.09, 34.75146), long = c(25.13, 29.45311, 14.0746, 14.45695, 
-2.916, 22.05, 74.44, 72.40194), date_BP = c(7429.833, 8048.077, 
4200, 4484.6, 4913.444, 8200.333, 3707.125, 2834.625)), row.names = c(NA, 
-8L), class = c("data.table", "data.frame"), .internal.selfref = <pointer: (nil)>)

Answer 2

Я хочу поблагодарить dc37 за помощь!

Я использовал stat_density2d, но получаю странные результаты.

Это первый подход:

world <- map_data("world")
plot <- ggplot()
plot <- ggplot(data=datgeo, aes(long,lat,fill=date_BP)) + 
  stat_density2d(aes(fill=..level..,alpha=..level..),geom='polygon',
                 colour='black') + 
  scale_fill_continuous(low="green",high="red") + guides(alpha="none")
plot <- plot + geom_map(dat=world, map = world, aes(map_id=region), 
                        fill="NA", color="black", inherit.aes = F)
plot <- plot + theme(panel.grid=element_blank(), 
                     panel.border=element_blank())
plot <- plot + coord_sf(xlim = c(-12.3, 110), ylim = c(70, 30), expand = FALSE)
windows()
plot

Я получил следующий сюжет:

Я не понимаю Масштаб легенды, так как по моим данным максимальное значение составляет 45020, а минимальное значение равно 736. Я не знаю, как это изменить.

Теперь, если я хочу добавить geom_point (), все становится странным:

world <- map_data("world")
plot <- ggplot()
plot <- ggplot(data=datgeo, aes(long,lat,fill=date_BP)) + 
  stat_density2d(aes(fill=..level..,alpha=..level..),geom='polygon',
                 colour='black') + 
  scale_fill_continuous(low="green",high="red") + guides(alpha="none")
plot <- plot + geom_map(dat=world, map = world, aes(map_id=region), 
                        fill="NA", color="black", inherit.aes = F)
plot <- plot + theme(panel.grid=element_blank(), 
                     panel.border=element_blank())
plot <- plot + coord_sf(xlim = c(-12.3, 110), ylim = c(70, 30), expand = FALSE)
plot <- plot + geom_point(data = datgeo, aes(x = long, y = lat))
windows()
plot

И я получаю эту карту:

Вы можете видеть, что теперь у меня есть точки, цвета изобар изменились и масштаб легенды в правильном формате. Я не понимаю, почему карта меняется таким образом только для добавления точек.

Знаете ли вы, что происходит, и какая карта является правильной?