Foxall's G-function with polygons in R spatstat

Question

I want to explore the distance of points in a point pattern to the nearest polygon in a polygon shapefile.

Reading the spatstat manual, vignettes and the book (chapter 8 in particular) I think I should to use the Gfox and Jfox functions in spatstat.

I imported the points and polygons shapefiles using package "sf" and converted to SpatialPointsDataframe, SpatialPolygonDataframe before converting the Points in a ppp object and the polygons in an owin object:

datafiles can be found here: https://www.dropbox.com/sh/ho8gp1rgpi0r7de/AABfvW3NdjinwjlZ0E5SRjGCa?dl=0

please adapt the code for importing the data according to where you store the files.

Load transect shapefile (the polygon containing observed points and polygons) and the polygons shapefiles

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

for Points I did Import dataframe and make it an sf object and convert it to ppp object

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)

and finally changed the Window to be sure it corresponded to the actual irregualr polygon from which the points were coming from (points sampled from a GIS layer as the intersection with a polygon layer)

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

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

I converted the polygons shapefile I needed to a Spatial object like:


rockT1_sp<-as(rockT1,'Spatial')

made it an owin object with:

rockT1_win<-as(rockT1_sp,'owin')

Finally tried to use Gfox:

fox<-Gfox(urc_poly,rockT1_win)

But I get the error:

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

Here some info on the objects:

> 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

I also thought to set the enclosing rectangle of rockT1_win equal to the one of urc_poly but I didn't find a way to to id.

Reading the example at page 283 in the book (section "8.10 Distance from a point pattern to another spatial pattern") it seem that should work. also the Gfox help says that Y should be an object of class "ppp", "psp" or "owin" to which distances will be measured.

Any idea of where there could be a problem?

PS: my 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

Thanks for a fairly detailed question. It would be much easier to diagnose this with a fully reproducible example. Can you provide that? — Ege Rubak
I would but i was limited in time... how can i provide you with the data, the shape file in particular? — Filippo
@EgeRubak I provided the data and modified the question so that the code should now work. (just adapt the data import part to your path). — Filippo

Ege Rubak Ege Rubak · Accepted Answer · 2019-08-09T23:05:01

Thanks for the fully reproducible example and an overall well written question.

I think you have stumbled into a bug in Gfox when Y is a owin object. I have made a quick fix in a branch on GitHub, but I need to think more about whether this is the best way to handle the problem, so the code may change in the near future.

If you want to try it out right away you will have to install from GitHub:

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

With this version of spatstat your code above runs without errors.

You write "I think I should to use the Gfox and Jfox functions in spatstat" and you may be right. It depends a bit on your scientific question and approach. If you simply want to investigate how the intensity of random points in X depends on the distance to the pattern Y when you condition on Y as being known and non-random, you could also use rhohat() and more formal techniques from Chapters 6 and 9 in the book.

Foxall's G-function with polygons in R spatstat

1 Answers