Foxall 的 G 函数与 R spatstat 中的多边形
Foxall's G-function with polygons in R spatstat
我想探索点模式中的点到多边形 shapefile 中最近的多边形的距离。
阅读 spatstat 手册、插图和书籍（特别是第 8 章）我认为我应该使用 spatstat 中的 Gfox 和 Jfox 函数。
我使用包 "sf" 导入了点和多边形 shapefile，并在转换 ppp 对象中的点和 owin 对象中的多边形之前转换为 SpatialPointsDataframe、SpatialPolygonDataframe：
数据文件可以在这里找到：https://www.dropbox.com/sh/ho8gp1rgpi0r7de/AABfvW3NdjinwjlZ0E5SRjGCa?dl=0
请根据您存储文件的位置修改导入数据的代码。
加载横断面形状文件（包含观察点和多边形的多边形）
和多边形 shapefile
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)
并最终更改 Window 以确保它对应于点来自的实际不规则多边形（从 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 转换为空间对象，例如：
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 应该是 class "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
感谢您提供完全可重现的示例和一个整体写得很好的问题。
我认为当 Y 是一个 owin 对象时，您在 Gfox 中遇到了一个错误。我已经在 GitHub 上的一个分支中进行了快速修复，但我需要更多地考虑这是否是处理问题的最佳方式，因此代码可能会在不久的将来发生变化。
如果您想立即试用，您必须从 GitHub:
安装
remotes::install_github("spatstat/spatstat", ref = "Gfox")
使用此版本的 spatstat 您上面的代码运行没有错误。
你写 "I think I should to use the Gfox and Jfox functions in spatstat" 你可能是对的。这在一定程度上取决于您的科学问题和方法。如果您只是想研究 X 中随机点的强度如何取决于与模式 Y 的距离，当您以 Y 为已知且非随机的条件时，您还可以使用 rhohat() 和章节中更正式的技术书中的6和9.
我想探索点模式中的点到多边形 shapefile 中最近的多边形的距离。
阅读 spatstat 手册、插图和书籍（特别是第 8 章）我认为我应该使用 spatstat 中的 Gfox 和 Jfox 函数。
我使用包 "sf" 导入了点和多边形 shapefile，并在转换 ppp 对象中的点和 owin 对象中的多边形之前转换为 SpatialPointsDataframe、SpatialPolygonDataframe：
数据文件可以在这里找到：https://www.dropbox.com/sh/ho8gp1rgpi0r7de/AABfvW3NdjinwjlZ0E5SRjGCa?dl=0
请根据您存储文件的位置修改导入数据的代码。
加载横断面形状文件（包含观察点和多边形的多边形） 和多边形 shapefile
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)
并最终更改 Window 以确保它对应于点来自的实际不规则多边形（从 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 转换为空间对象，例如：
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 应该是 class "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
感谢您提供完全可重现的示例和一个整体写得很好的问题。
我认为当 Y 是一个 owin 对象时，您在 Gfox 中遇到了一个错误。我已经在 GitHub 上的一个分支中进行了快速修复，但我需要更多地考虑这是否是处理问题的最佳方式，因此代码可能会在不久的将来发生变化。
如果您想立即试用，您必须从 GitHub:
安装remotes::install_github("spatstat/spatstat", ref = "Gfox")
使用此版本的 spatstat 您上面的代码运行没有错误。
你写 "I think I should to use the Gfox and Jfox functions in spatstat" 你可能是对的。这在一定程度上取决于您的科学问题和方法。如果您只是想研究 X 中随机点的强度如何取决于与模式 Y 的距离，当您以 Y 为已知且非随机的条件时，您还可以使用 rhohat() 和章节中更正式的技术书中的6和9.