R 用空间数据的值填充矩阵

Question

我在 R 中有一个名为“habitat”的矩阵，其中纬度作为行名，经度作为列名。细胞现在充满了 NA。我的最终目标是让整个矩阵填充值 - 0 或来自另一个数据帧的值。

对于陆地上的单元格，我希望单元格值为 0。

我使用 R 中的 sf 包创建了一个简单的要素集合，其中的边界框与栖息地矩阵的纬度和经度相匹配，以及 1 个土地多边形要素。

land<-st_read(paste0(nc_data_path, land))
land <- st_transform(land, crs = "EPSG:4326")
vertices <- matrix(c(bbox["xmin"], bbox["ymin"], 
                     bbox["xmax"], bbox["ymin"], 
                     bbox["xmax"], bbox["ymax"], 
                     bbox["xmin"], bbox["ymax"], 
                     bbox["xmin"], bbox["ymin"]), ncol = 2, byrow = TRUE)
polygon <- st_polygon(list(vertices))
polygon <- st_sfc(polygon)
st_crs(polygon) <- st_crs(st_crs("EPSG:4326"))
land=st_crop(land, polygon) 

我不知道如何使用土地简单特征在有土地的地方用0填充栖息地矩阵。我尝试将栖息地矩阵转换为空间点数据框，但我不确定从哪里开始

latitudes <- as.numeric(rownames(habitat))
longitudes <- as.numeric(colnames(habitat))
habitat_dt <- data.table(latitude = rep(latitudes, each = ncol(habitat)),
                         longitude = rep(longitudes, times = nrow(habitat)),
                         value = as.vector(habitat))
habitat_points_sf <- st_as_sf(habitat_dt, coords = c("longitude", "latitude"), crs = 4326)

对于不在陆地上（在海洋中）的单元格，我希望根据最接近的纬度/经度（sdm 数据框列是纬度、经度和磅）从名为“sdm”的数据帧填充值。我想避免使用 for 循环来执行此操作。

Answer 1

解决方案和说明如下（其中

habitats_dt

是带有几何列的 data.table，

sdm_sf

是

land_sf

是 sf 数据框）。

使用

sf::st_within

标记海洋栖息地：

habitats_dt[, ocean := is.na(as.logical(st_within(geometry, land_sf)))]

使用

nngeo::st_nn

:

将 sdm 最近邻居的行 ID 添加到每个海洋栖息地

habitats_dt[ocean==TRUE, sdm_rowid := transpose(st_nn(geometry, sdm_sf))]

从 sdm 数据添加相应的 lbs 属性（0 表示陆地单元）

habitats_dt[, lbs := sdm_sf$lbs[c(habitats_dt$sdm_rowid)]][is.na(lbs), lbs := 0]

转换为矩阵：

result <- as.matrix(dcast(habitats_dt, lat ~ lon, value.var="lbs"), rownames=T)

图解：

library(data.table)
library(sf)
library(nngeo)

# DUMMY DATA
# habitats
min_lat <- -4; max_lat <- min_lat + 10; min_lon <- 40; max_lon <- min_lon + 10
lat <- seq(min_lat, max_lat, 1) # as.numeric(rownames(habitat))
lon <- seq(min_lon, max_lon, 2) # as.numeric(colnames(habitat))
habitats_dt <- data.frame(lat = rep(lat, each = length(lon)), lon = rep(lon, times = length(lat))) |>
  st_as_sf(coords = c("lon", "lat"), crs = 4326, remove=F) |>
  setDT()
# land
land_sf <- st_polygon(list(cbind(min_lon+c(1.5,3.5,9.5,1.5), min_lat+c(1.5,9.5,6.5,1.5))))
land_sf <- st_sf(geometry = st_sfc(land_sf), crs = 4326)
# sdm
set.seed(123)
n <- 15
sdm_sf <- data.frame(
  lat = runif(n, min = min_lat, max = max_lat),
  lon = runif(n, min = min_lon, max = max_lon),
  lbs = seq(100,999,length.out=n)) |>
  st_as_sf(coords = c("lon", "lat"), crs = 4326)
sdm_sf <- setDT(sdm_sf)[is.na(as.logical(st_within(geometry, land_sf)))] |> st_as_sf()

class(habitats_dt)
class(land_sf)
class(sdm_sf)

# SOLUTION
habitats_dt[, ocean := is.na(as.logical(st_within(geometry, land_sf)))]
habitats_dt[ocean==TRUE, sdm_rowid := transpose(st_nn(geometry, sdm_sf))]
habitats_dt[, lbs := sdm_sf$lbs[c(habitats_dt$sdm_rowid)]][ocean==FALSE, lbs := 0]
result <- as.matrix(dcast(habitats_dt, lat ~ lon, value.var="lbs"), rownames=T)   

# VISUAL CHECK
print(result)
library(ggplot2)
ggplot() + 
  geom_sf(data=land_sf, fill="grey") +
  geom_sf(data=st_as_sf(habitats_dt), aes(colour=lbs, size=lbs)) +
  geom_sf(data=sdm_sf, aes(colour=lbs, size=lbs), shape="X") +
  guides(size="none")