使用栅格网格提取行数据

Question

我有一个分辨率为 0.5 度 (r) 的栅格网格和一个包含 3 列的数据框 (my_df)：long、lat 和 id。数据框代表物种出现记录。

我想要做的是确定我的栅格网格的每个 0.5 度单元格中存在哪些物种，并且每个单元格只保留每个物种的 1 条记录（my_df 有超过 90,000,000 行），所以如果一个 0.5 度单元格只有一个物种，会有一行包含栅格网格单元的纬度、经度，然后是数据框中的物种 ID。其他栅格网格单元可能包含数百个物种，因此可能有数百行。

最终我想创建一个数据框，其中包含每个物种位置所在的 0.5 度栅格网格的长和纬度以及那里存在的物种 ID，每个物种一行。

我已经创建了一个栅格网格，按照...

ext <- extent(-180.0, 180, -90.0, 90.0)
gridsize <- 0.5
r <- raster(ext, res=gridsize)
crs(r) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

还有一个数据框，它原本是一个 SpatialPolygonsDataframe...

A tibble: 6 x 3
   long   lat id   
  <dbl> <dbl> <chr>
1  16.5 -28.6 0    
2  16.5 -28.6 0    
3  16.5 -28.6 0    
4  16.5 -28.6 0    
5  16.5 -28.6 0    
6  16.5 -28.6 0 
etc
etc

...但我不确定如何继续该方法的其余部分。我曾尝试栅格化我的数据、提取点等，但我不断遇到错误并且不确定用于实现我的目标的正确方法。

或者，如果有人知道如何直接从 SpatialPolygonsDataFrame 中提取物种名称，其中包含每个物种的范围多边形，位于 0.5 度栅格网格单元位置，那就太好了。

任何帮助将不胜感激。

Answer 1

使用点数据，您可以这样做

示例数据

#species
set.seed(0)
n <- 20
spp <- data.frame(lon=runif(n, -180, 180), lat=runif(n,-90,90), sp=sample(5, n, replace=TRUE)) 

# raster
library(raster)
# for the example I use a resolution of 90, rather than 0.5 
r <- raster(res=90)

现在计算每个位置的单元格编号并制成表格。我这样做的方式是返回计数，而不仅仅是存在/不存在

spp$cell <- cellFromXY(r, spp[, c("lon", "lat")])
tb <- table(spp$cell, spp$sp)

获取每个单元格的 lon/lat

xy <- xyFromCell(r, as.integer(rownames(tb)))
result <- cbind(xy, tb)
colnames(result)[1:2] <- c("lon", "lat")
result
#   lon lat 1 2 3 4 5
#1 -135  45 0 0 1 0 0
#2  -45  45 0 2 1 0 0
#3   45  45 1 0 0 2 0
#4  135  45 0 1 0 0 1
#5 -135 -45 1 2 0 0 0
#6  -45 -45 0 1 0 1 0
#7   45 -45 1 1 0 0 0
#8  135 -45 1 0 1 2 0

对于多边形数据（以及点数据），您可以使用raster::rasterize

多边形数据示例

library(raster)
p1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))
hole <- rbind(c(-150,-20), c(-100,-10), c(-110,20), c(-150,-20))
p1 <- list(p1, hole)
p2 <- rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0))
p3 <- rbind(c(-125,0), c(0,60), c(40,5), c(15,-45), c(-125,0))
spp <- data.frame(species=letters[1:3], stringsAsFactors=FALSE)
pols <- spPolygons(p1, p2, p3, attr=spp)

对每个物种进行栅格化并组合在一个 RasterStack 中。如果您有许多物种，您想为 rasterize 参数分配一个文件名，例如 filename = paste0("sp_", i, ".tif")

usp <- unique(spp$species)
r <- raster(res=0.5)
s <- list()
for (i in 1:length(usp)) {
    p <- pols[pols$species == usp[i], ]
    s[[i]] <- rasterize(p, r, field=1, fun="count")
}       
ss <- stack(s)

（对于物种丰富度做sr <- sum(ss>0, na.rm=TRUE)）

创建你想要的输出

m <- as.matrix(ss)
m[is.na(m)] <- 0
# to remove rows with no species 
i <- which(rowSums(m) > 0)
xy <- xyFromCell(r, i)  
output <- cbind(xy, m[i,])
colnames(output) <- c("lon", "lat", usp)
head(output)
#        lon   lat a b c
#[1,]  -0.25 59.75 0 0 1
#[2,] 139.75 59.75 0 1 0
#[3,]  -1.25 59.25 0 0 1
#[4,]  -0.75 59.25 0 0 1
#[5,]  -0.25 59.25 0 0 1
#[6,]   0.25 59.25 0 0 1

Answer 2

如果我猜对了，您想匹配单元格内的点。我认为您正在寻找基于点和多边形之间相交的空间连接。

我强烈建议您使用sf 包而不是sp 对象。这就是我要向你提出的建议。

首先，使用st_make_grid函数创建网格

library(sf)
library(dplyr)

ext <- raster::extent(-180.0, 180, -90.0, 90.0)

grid <- st_bbox(ext) %>% 
  st_make_grid(cellsize = 0.5, what = "polygons") %>%
  st_set_crs(4326)
grid <- grid %>% st_sf() %>% mutate(id_cell = seq_len(nrow(.)))

那么我们来做一个简单的dataframe：

df <- data.frame(long = 16.51, lat = -28.6, id = 0)
df <- df %>% sf::st_as_sf(coords = c("long","lat"), crs = 4326)

df

Simple feature collection with 1 feature and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: 16.51 ymin: -28.6 xmax: 16.51 ymax: -28.6
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
  id            geometry
1  0 POINT (16.51 -28.6)

然后，您需要使用st_join 函数。默认情况下，空间连接基于交集：

df %>% sf::st_join(grid, left = TRUE)

although coordinates are longitude/latitude, st_intersects assumes that they are planar
Simple feature collection with 1 feature and 2 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: 16.51 ymin: -28.6 xmax: 16.51 ymax: -28.6
epsg (SRID):    4326
proj4string:    +proj=longlat +datum=WGS84 +no_defs
  id id_cell            geometry
1  0   88234 POINT (16.51 -28.6)

我假设你想要一个左连接（报告你的所有观点）。您可以更改该选项。我认为使用sf 会比手动编码技术更快。