r - 创建密度栅格并按多边形特征提取总和

Question

我有一个多边形( zones )和一组坐标( points )。我想为整个多边形创建一个空间核密度栅格，并按区域提取密度总和。多边形之外的点应该被丢弃。

library(raster)
library(tidyverse)
library(sf)
library(spatstat)
library(maptools)

load(url("https://www.dropbox.com/s/iv1s5butsx2v01r/example.RData?dl=1"))
# alternatively, links to gists for each object
# https://gist.github.com/ericpgreen/d80665d22dfa1c05607e75b8d2163b84
# https://gist.github.com/ericpgreen/7f4d3cee3eb5efed5486f7f713306e96

ggplot() +
  geom_sf(data = zones) +
  geom_sf(data = points) +
  theme_minimal()

我尝试转换为 ppp与 { spatstat然后使用 density() ，但我对结果中的单位感到困惑。我相信问题与 map 的单位有关，但我不确定如何进行。

更新

这是重现我创建的密度图的代码:

zones_owin <- as.owin(as_Spatial(zones))
pts <- st_coordinates(points)
p <- ppp(pts[,1], pts[,2], window=zones_owin, unitname=c("metre","metres"))
ds <- density(p) 
r <- raster(ds)
plot(r)

Answer 1

当您直接使用地理坐标(经度、纬度)时，单位是困难的。如果可能，您应该转换为平面坐标(这是 spatstat 的要求)并从那里继续。平面坐标通常以米为单位，但我想这取决于具体的投影和底层椭球等。你可以看到 this answer关于如何使用 sf 投影到平面坐标并导出到 spatstat格式使用 maptools . 注意:您必须手动选择一个合理的投影(您可以使用 http://epsg.io 找到一个)并且您必须同时投影多边形和点。

一旦一切都在 spatstat您可以使用的格式 density.ppp做核平滑。生成的网格值(类 im 的对象)是点的强度，即每平方单位(例如平方米)的点数。如果你想聚合某个区域，你可以使用 integral.im(..., domain = ...)获得具有给定强度的点过程模型在该区域中的预期点数。