如何从“ppm”的地理栅格中获取协变量数据?
How to get covariate data from a geographic raster for `ppm`?
我想用 spatstat::ppm
拟合泊松点过程模型,但我不确定将协变量数据提供给函数的最佳方法是什么。我知道 spatstat
需要平面坐标,所以我在创建 ppp
点模式对象之前将我的点位置数据转换为平面 crs。协变量数据位于具有未投影地理坐标的栅格堆栈中,据我所知 projecting rasters is generally ill-advised。我使用点的原始地理坐标和 raster::extract
从栅格中提取点位置的协变量值。到目前为止,一切都很好。问题是...
it is not sufficient to have observed the covariate only at the points
of the data point pattern; the covariate must also have been observed
at other locations in the window. -ppm
helpfile
我似乎有两个选项可以为 data
参数提供协变量数据。
- 一张像素图片;由于光栅投影问题,这似乎是不明智的。
- 可以在任何位置 (x,y) 进行评估以获得相应的协变量值的函数列表(每个协变量一个)。这似乎是要走的路,但我尝试编写这样一个函数的速度却慢得离谱。在将坐标转换为栅格的 crs 之后,它会为每个坐标对调用
raster::extract
。虽然 raster::extract
在给定大量点时速度相当快,但每次调用似乎都会产生大量开销。根据microbenchmark
,坐标变换大约需要 4ms,提取单个协变量大约需要 582ms,或者每个点大约需要 4 秒才能得到所有 7 个协变量。我不知道 ppm
会调用多少次,但如果模式中的每个点调用一次,它会花费太长时间。
有没有什么方法可以找出 ppm
将查询协变量数据的完整点集,以便我可以通过一次调用预先提取这些点?
似乎我的用例(地理栅格中的协变量)应该很常见,所以我猜有一种既定的方法可以正确地做到这一点。这是什么?
感谢您提出一个写得很好的问题,清楚地表明您的需求。使用例如一个简单的可重现示例会更好来自 raster
和 spatstat
的内置数据或人工生成的数据。由于缺少可重现的示例,我的答案将不包含任何代码,但会概述您可以做什么。
ppm
中的第一步是制定正交方案或 class quad
或 logiquad
,具体取决于 ppm
中使用的最大似然逼近.这些可以由用户通过 quadscheme
或 quadscheme.logi
直接生成。正交方案包含 ppm
将评估协变量的所有点。您可以使用函数 coords
提取正交方案的坐标。如果您构建一个 data.frame
并在这些点评估所有协变量,您可以将其作为 data
参数提供给 ppm
,而正交方案是第一个参数。要更好地理解事物,请尝试阅读 help(ppm.quad)
.
的详细信息部分
另一种可以优化数据使用的方法是提取当前栅格堆栈的网格点以及所有协变量值并投影该点数据。然后将其转换为简单的 data.frame
,其中包含列 x
、y
、covar1
、covar2
等。然后您可以使用 x
和 y
连同您感兴趣的观察点以手动创建正交方案,其余列可以作为 data
到 ppm
.
提供
比较这两种方法的结果以及仅投影栅格堆栈并将其转换为 im
对象列表的结果会很有趣。
我想用 spatstat::ppm
拟合泊松点过程模型,但我不确定将协变量数据提供给函数的最佳方法是什么。我知道 spatstat
需要平面坐标,所以我在创建 ppp
点模式对象之前将我的点位置数据转换为平面 crs。协变量数据位于具有未投影地理坐标的栅格堆栈中,据我所知 projecting rasters is generally ill-advised。我使用点的原始地理坐标和 raster::extract
从栅格中提取点位置的协变量值。到目前为止,一切都很好。问题是...
it is not sufficient to have observed the covariate only at the points of the data point pattern; the covariate must also have been observed at other locations in the window. -
ppm
helpfile
我似乎有两个选项可以为 data
参数提供协变量数据。
- 一张像素图片;由于光栅投影问题,这似乎是不明智的。
- 可以在任何位置 (x,y) 进行评估以获得相应的协变量值的函数列表(每个协变量一个)。这似乎是要走的路,但我尝试编写这样一个函数的速度却慢得离谱。在将坐标转换为栅格的 crs 之后,它会为每个坐标对调用
raster::extract
。虽然raster::extract
在给定大量点时速度相当快,但每次调用似乎都会产生大量开销。根据microbenchmark
,坐标变换大约需要 4ms,提取单个协变量大约需要 582ms,或者每个点大约需要 4 秒才能得到所有 7 个协变量。我不知道ppm
会调用多少次,但如果模式中的每个点调用一次,它会花费太长时间。
有没有什么方法可以找出 ppm
将查询协变量数据的完整点集,以便我可以通过一次调用预先提取这些点?
似乎我的用例(地理栅格中的协变量)应该很常见,所以我猜有一种既定的方法可以正确地做到这一点。这是什么?
感谢您提出一个写得很好的问题,清楚地表明您的需求。使用例如一个简单的可重现示例会更好来自 raster
和 spatstat
的内置数据或人工生成的数据。由于缺少可重现的示例,我的答案将不包含任何代码,但会概述您可以做什么。
ppm
中的第一步是制定正交方案或 class quad
或 logiquad
,具体取决于 ppm
中使用的最大似然逼近.这些可以由用户通过 quadscheme
或 quadscheme.logi
直接生成。正交方案包含 ppm
将评估协变量的所有点。您可以使用函数 coords
提取正交方案的坐标。如果您构建一个 data.frame
并在这些点评估所有协变量,您可以将其作为 data
参数提供给 ppm
,而正交方案是第一个参数。要更好地理解事物,请尝试阅读 help(ppm.quad)
.
另一种可以优化数据使用的方法是提取当前栅格堆栈的网格点以及所有协变量值并投影该点数据。然后将其转换为简单的 data.frame
,其中包含列 x
、y
、covar1
、covar2
等。然后您可以使用 x
和 y
连同您感兴趣的观察点以手动创建正交方案,其余列可以作为 data
到 ppm
.
比较这两种方法的结果以及仅投影栅格堆栈并将其转换为 im
对象列表的结果会很有趣。