从空间数据制作网格图

Question

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

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

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

Answer 1

这是我使用 sf 包的尝试。首先，我导入了您的数据并将其转换为 sf 对象。然后，我创建了另一个包含网格的 sf 对象。我使用 raster 包和 sf 包来创建网格。有了两个 sf 对象后，我计算了每个网格中存在多少数据点，并将结果作为新列添加到 foo 中。最后，我画了一个图形。

library(tidyverse)
library(sf)
library(raster)
library(viridis)

# Import the data and convert it to an sf object
mydata <- read_csv("https://raw.githubusercontent.com/abenmhamed/data/main/Events.csv") %>%
          st_as_sf(coords = c("LON", "LAT"), 
                   crs = 4326, agr = "constant")

# Create an sf object for the grid
gridx <- seq(from = -19,to = -10, by = 0.083333)
gridy <- seq(from = 20,to = 29, by = 0.083333)


foo <- raster(xmn = -19, xmx = -10,
              ymn = 20, ymx = 29,
              nrows = length(gridx),
              ncols = length(gridy)) %>% 
       rasterToPolygons() %>% 
       st_as_sf(crs = 4326) %>% 
       mutate(group = 1:(length(gridx)*length(gridy))) %>% 
       st_cast("MULTIPOLYGON")


# Now count how many data points exist in each grid
mutate(foo,
       count = lengths(st_intersects(x = foo,  y = mydata))) -> foo

# Draw a graphic
ggplot() +
  geom_sf(data = foo, aes(fill = count)) +
  scale_fill_viridis(option = "D") -> g

Answer 2

不错的数据集，假设这些是渔船 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))