Foxall 的 G 函数在 R spatstat 中具有多边形

Posted 2023-03-23

【中文标题】Foxall 的 G 函数在 R spatstat 中具有多边形

【英文标题】：Foxall's G-function with polygons in R spatstat

【发布时间】：2019-12-16 16:39:03

【问题描述】：

我想探索点模式中的点到多边形 shapefile 中最近的多边形的距离。

阅读 spatstat 手册、小插图和书（特别是第 8 章）我认为我应该在 spatstat 中使用 Gfox 和 Jfox 函数。

在转换 ppp 对象中的点和 owin 对象中的多边形之前，我使用包“sf”导入了点和多边形 shapefile 并转换为 SpatialPointsDataframe、SpatialPolygonDataframe：

数据文件可以在这里找到：https://www.dropbox.com/sh/ho8gp1rgpi0r7de/AABfvW3NdjinwjlZ0E5SRjGCa?dl=0

请根据您存储文件的位置调整用于导入数据的代码。

加载样带形状文件（包含观察点和多边形的多边形） 和多边形形状文件

library(tidyverse)
 # load these two together because spatstat rely on them but I don't know exactly for what.
library(sp)
library(maptools)  # needed for method such as as.ppp

library(spatstat)
library(sf)

trs    <-st_read('transect.shp')
rockT1 <-st_read('polygons.shp')  # shapefiles containing the polygons from which I want to calculate the distance to the points using Gfox

对于我所做的点 导入dataframe并使其成为sf对象并转换为ppp对象

points<-read.table('points.txt', head=T,sep='\t',dec='.') 

options(digits=15) # this to allow enough decimals
urcT1_sp<-points %>% 
        mutate(geometry=as.character(geometry)) %>%
        mutate(lon=as.numeric(sapply(strsplit(geometry, '[(,)]'), "[[", 2)), 
               lat=as.numeric(sapply(strsplit(geometry, '[(,)]'), "[[", 3))) %>%
        st_as_sf(coords=c('lon','lat'),crs= 32619) %>% 
        as(.,'Spatial')

urc_poly<-as.ppp(urcT1_sp)

最后更改了窗口以确保它对应于点来自的实际不规则多边形（从 GIS 层采样的点作为与多边形层的交集）

tr_poly_sp<- trs %>% select(geometry) %>% as(., 'Spatial') %>% as(.,'SpatialPolygons')

tr_poly_win<- as.owin(tr_poly_sp)
Window(urc_poly)<-tr_poly_win

我将需要的多边形 shapefile 转换为 Spatial 对象，例如：

rockT1_sp<-as(rockT1,'Spatial')

通过以下方式使其成为 owin 对象：

rockT1_win<-as(rockT1_sp,'owin')

终于尝试用Gfox了：

fox<-Gfox(urc_poly,rockT1_win)

但我得到了错误：

Error in owin(range(x), range(y)) : 
 xrange should be a vector of length 2 giving (xmin, xmax)
In addition: Warning messages:
1: In Gfox(urc_poly, rockT1_win) :
 Trimming the window of X to be a subset of the window of Y
2: In min(x) : no non-missing arguments to min; returning Inf
3: In max(x) : no non-missing arguments to max; returning -Inf

这里有一些关于对象的信息：

> urc_poly
Marked planar point pattern: 98 points
Mark variables: unique_ID, site, transect, pos, sub_type, id_merged, study 
window: polygonal boundary
enclosing rectangle: [687281.931393065, 687309.776978588] x [5559021.04635051, 5559034.27763186] units

> rockT1_win
window: polygonal boundary
enclosing rectangle: [687286.164022728, 687309.715980057] x [5559023.59638782, 5559033.57055627] units

我还想将 rockT1_win 的封闭矩形设置为 urc_poly 的一个，但我没有找到 id 的方法。

阅读本书第 283 页的示例（“8.10 点模式到另一个空间模式的距离”部分）似乎应该可行。 Gfox 的帮助还说 Y 应该是“ppp”、“psp”或“owin”类的对象，距离将被测量。

知道哪里有问题吗？

PS：我的 sessionInfo()

R version 3.6.0 (2019-04-26)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 17763)

Matrix products: default

locale:
[1] LC_COLLATE=Italian_Italy.1252  LC_CTYPE=Italian_Italy.1252    LC_MONETARY=Italian_Italy.1252 LC_NUMERIC=C                   LC_TIME=Italian_Italy.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] xlsx_0.6.1          sf_0.7-4            spatstat_1.60-1     rpart_4.1-15        nlme_3.1-140        spatstat.data_1.4-0 maptools_0.9-5      sp_1.3-1            forcats_0.4.0       stringr_1.4.0       dplyr_0.8.1        
[12] purrr_0.3.2         readr_1.3.1         tidyr_0.8.3         tibble_2.1.1        ggplot2_3.1.1       tidyverse_1.2.1    

loaded via a namespace (and not attached):
[1] Rcpp_1.0.1            lubridate_1.7.4       lattice_0.20-38       deldir_0.1-16         xlsxjars_0.6.1        class_7.3-15          digest_0.6.18         assertthat_0.2.1      R6_2.4.0              cellranger_1.1.0     
[11] plyr_1.8.4            backports_1.1.4       evaluate_0.13         e1071_1.7-1           httr_1.4.0            tensor_1.5            pillar_1.4.0          rlang_0.3.4           lazyeval_0.2.2        readxl_1.3.1         
[21] rstudioapi_0.10       Matrix_1.2-17         goftest_1.1-1         rmarkdown_1.12        splines_3.6.0         rgdal_1.4-3           foreign_0.8-71        polyclip_1.10-0       munsell_0.5.0         broom_0.5.2          
[31] xfun_0.7              compiler_3.6.0        modelr_0.1.4          pkgconfig_2.0.2       mgcv_1.8-28           htmltools_0.3.6       tidyselect_0.2.5      crayon_1.3.4          withr_2.1.2           grid_3.6.0           
[41] jsonlite_1.6          gtable_0.3.0          DBI_1.0.0             magrittr_1.5          units_0.6-3           scales_1.0.0          KernSmooth_2.23-15    cli_1.1.0             stringi_1.4.3         xml2_1.2.0           
[51] spatstat.utils_1.13-0 generics_0.0.2        tools_3.6.0           glue_1.3.1            hms_0.4.2             abind_1.4-5           colorspace_1.4-1      classInt_0.3-3        rvest_0.3.4           rJava_0.9-11         
[61] knitr_1.22            haven_2.1.0 

【问题讨论】：

感谢您提出相当详细的问题。使用完全可重现的示例来诊断这一点会容易得多。你能提供吗？

我会，但我时间有限......我如何为您提供数据，特别是形状文件？

@EgeRubak 我提供了数据并修改了问题，以便代码现在可以工作。 （只需根据您的路径调整数据导入部分）。

@EgeRubak 刚刚测试了代码，错误重现。

【参考方案1】：

感谢您提供完全可重复的示例和一个整体写得很好的问题。

当Y 是owin 对象时，我认为您偶然发现了Gfox 中的错误。我已经在 GitHub 上的一个分支中进行了快速修复，但我需要更多地考虑这是否是处理问题的最佳方式，因此代码可能会在不久的将来发生变化。

如果您想立即试用，则必须从 GitHub 安装：

remotes::install_github("spatstat/spatstat", ref = "Gfox")

使用此版本的spatstat，您的上述代码可以正常运行。

您写道“我认为我应该在 spatstat 中使用 Gfox 和 Jfox 函数”，您可能是对的。这在一定程度上取决于您的科学问题和方法。如果您只是想研究当您将 Y 设为已知且非随机时，X 中随机点的强度如何取决于与模式 Y 的距离，您还可以使用rhohat() 和第 6 章中的更正式的技术和9 本书。

【讨论】：

非常感谢您的宝贵时间和回答！

我对空间统计很陌生，所以我正在探索我能做什么。事实上，在阅读了 Fortin 和 Dale 2002 上的 G 距离后，我直接进入了第 8 章。我很快就会阅读更多的书和第 6 章和第 9 章！感谢您的建议！