Извлечение данных строки с использованием растровой сетки - PullRequest
Новая
       71

Извлечение данных строки с использованием растровой сетки

0 голосов
/ 07 апреля 2020

У меня есть растровая сетка с разрешением 0,5 градуса (r) и фрейм данных (my_df) с 3 столбцами: long, lat и id. Фрейм данных представляет записи встречаемости видов.

Что я хочу сделать, это определить, какие виды присутствуют в каждой ячейке 0,5 градуса моей растровой сетки, и для каждой ячейки вести только 1 запись каждого вида (my_df имеет более 90 000 000 строк), поэтому, если ячейка 0,5 градуса имеет только один вид, будет строка с широтой, длиной ячейки растровой сетки и идентификатором вида из кадра данных. Другие ячейки растровой сетки могут содержать сотни видов, поэтому могут иметь сотни строк.

В конечном итоге я хотел бы создать фрейм данных, который имеет длину и широту растровой сетки 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

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

В качестве альтернативы, если кто-то знает, как извлечь названия видов напрямую сформируйте SpatialPolygonsDataFrame, который содержит полигон диапазона для каждого вида, в ячейках растровой сетки на 0,5 градуса, что было бы превосходно.

Любая помощь будет принята с благодарностью.

Ответы [ 2 ]

1 голос
/ 07 апреля 2020

Если я правильно угадал, вы хотите сопоставить точки, попадающие в ячейки. Я думаю, что вы ищете пространственное соединение, основанное на пересечении точек и многоугольников.

Я настоятельно рекомендую вам использовать пакет 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(.)))

Затем давайте возьмем простой кадр данных: 

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 будет быстрее, чем ручная кодировка.

0 голосов
/ 07 апреля 2020

С точечными данными вы можете сделать это следующим образом

Пример данных 

#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. Если у вас много видов, вы хотите присвоить имя файла аргументу растеризации, например 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
...