r - 从空间数据制作网格图

Question

我在数据框中有空间坐标，其中每一行（经度、纬度）对应于我正在关注的事件的发生。我试图映射这些数据，但我不想使用点，而是想创建一个网格，其单元格的分辨率为 5 海里（~ 0.083333），并计算每个单元格中事件的发生次数并绘制它。

这是我在一些资源的帮助下编写的代码。但它看起来不像我预期的那样。你能弄清楚我做错了什么吗？我附上了原始位置和我得到的结果地图。这是数据的链接。

re_pi = read.csv(file = "~/Desktop/Events.csv")

gridx <- seq(from=-19,to=-10,by=0.083333)
gridy <- seq(from=20,to=29,by=0.083333)
xcell <- unlist(lapply(re_pi$LON,function(x) min(which(gridx>x))))
ycell <- unlist(lapply(re_pi$LAT,function(y) min(which(gridy>y))))
re_pi$cell <- (length(gridx) - 1) * ycell + xcell

rr = re_pi %>%
  group_by(cell)%>%
  summarise(Lat = mean(LAT),Lon = mean(LON),Freq = length(cell))

my_theme <- theme_bw() + theme(panel.ontop=TRUE, panel.background=element_blank())
my_cols <- scale_color_distiller(palette='Spectral')
my_fill <- scale_fill_distiller(palette='Spectral')

ggplot(rr, aes(y=Lat, x=Lon, fill=Effort)) + geom_tile(width=1.2, height=1.2) +
  borders('world', xlim=range(rr$Lon), ylim=range(rr$Lat), colour='black') + my_theme + my_fill +
  coord_quickmap(xlim=range(rr$Lon), ylim=range(rr$Lat)) 

标签： rmapsspatialsf

Answer 1

不错的数据集，假设这些是渔船 VMS 数据。这可能是实现您的目标的一种方法，它严重依赖于 tidyverse 和绕过光栅和形状。

library(tidyverse)
library(mapdata) # higher resolution maps
# poor man's gridding function
grade <- function (x, dx) {
  if (dx > 1)
    warning("Not tested for grids larger than one")
  brks <- seq(floor(min(x)), ceiling(max(x)), dx)
  ints <- findInterval(x, brks, all.inside = TRUE)
  x <- (brks[ints] + brks[ints + 1])/2
  return(x)
}
d <-
  read_csv("https://raw.githubusercontent.com/abenmhamed/data/main/Events.csv") %>%
  janitor::clean_names() %>%
  # make a grid 0.01 x 0.01 longitude / latitude
  mutate(lon = grade(lon, 0.01),
         lat = grade(lat, 0.01)) %>%
  group_by(lon, lat) %>%
  count() %>%
  # not much happening south of 21 and north of 26
  filter(between(lat, 21, 26.25))
d %>%
  ggplot() +
  theme_bw() +
  geom_tile(aes(lon, lat, fill = n)) +
  scale_fill_viridis_c(option = "B", direction = -1) +
  # only data within the data-bounds
  borders(database = "worldHires",
          xlim = range(d$lon), ylim = range(d$lat),
          fill = "grey") +
  labs(x = NULL, y = NULL, fill = "Effort") +
  # limit plot
  coord_quickmap(xlim = range(d$lon), ylim = range(d$lat)) +
  # legends within plot
  theme(legend.position = c(0.77, 0.26))