Use sf polygon object as window in spatstat

Question

Hello all potential helpers,

I have a SpatialPolygonDataFrame object obtained from the tigris package and I would like to use it as a polygonal window in the creation of a ppp object. Here is what I tried:

# Needed packages
library(spatstat)
library(sf)

# Download geospatial data for Lee county in Alabama (home of the great Auburn University by the way!)
county <- tigris::county_subdivisions(state = "Alabama", county = "Lee")

# The polygon of Lee county is subdivided, so I convert it to a single polygon after converting it to an sf object
county_sf <- st_as_sf(county)
county_one <- st_union(county_sf)

# A quick plot of the object outputs what I am expecting
plot(county_one)

# Now I create a planar point pattern and I use county_one as the window
p <- ppp(x = -85.4, y = 32.5, window = as.owin((county_one)))

# But the plot here shows that the window is just a rectangle and not the polygon :(
plot(p)

Thank you for your help.

This is just in the making. See https://github.com/r-spatial/sf/issues/1233 for more information — TimSalabim, Jan 05 '20 at 07:27

Ege Rubak · Accepted Answer · 2020-01-06T16:50:05.833

Note: I have edited this answer to contain full details.

As @TimSalabim mentions this is under way in sf, but until then you have to go through the old sp classes such as SpatialPolygons. Use something like as_Spatial in sf and then load maptools and use as.owin or as(x, "owin") on the Spatial object.

Furthermore, you can only use coordinates in planar (projected) space with spatstat and not coordinates on the curved surface of the earth. You have to project to a relevant planar coordinates system. Maybe <epsg.io/6345> is usable in this case. To project to this coordinate system use sf::st_transform(county_one, crs = 6345). Afterwards you convert to Spatial and then owin. Note: Choosing the relevant projection is a science, and I don’t know much about it, so do a bit of research if you want to make sure you don’t get too distorted results.

Specifically with the original example you can do:

# Needed packages
library(spatstat)
#> Loading required package: spatstat.data
#> Loading required package: nlme
#> Loading required package: rpart
#> 
#> spatstat 1.62-2       (nickname: 'Shape-shifting lizard') 
#> For an introduction to spatstat, type 'beginner'
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.2, PROJ 6.2.1
library(maptools)
#> Loading required package: sp
#> Checking rgeos availability: TRUE
library(tigris)
#> To enable 
#> caching of data, set `options(tigris_use_cache = TRUE)` in your R script or .Rprofile.
#> 
#> Attaching package: 'tigris'
#> The following object is masked from 'package:graphics':
#> 
#>     plot

county <- county_subdivisions(state = "Alabama", county = "Lee", class = "sf", progress_bar = FALSE)
county_one <- st_union(county)
plot(county_one)

county_flat <- st_transform(county_one, crs = 6345)
plot(county_flat)

county_owin <- as.owin(as_Spatial(county_flat))

100 random points in the county:

p <- runifpoint(100, win = county_owin)
plot(p)

Thank you for your answer. The `sf` function is `as_Spatial`. I unfortunately get this error when I use `as.owin`: `Error in maptools:::as.owin.SpatialPolygons(county_sp) : Only projected coordinates may be converted to spatstat class objects`. — SavedByJESUS, Jan 05 '20 at 21:36
You can only use coordinates in planar (projected) space with `spatstat` and not coordinates on the curved surface of the earth. You have to project to relevant coordinates system. Maybe https://epsg.io/6345 is relevant. To project to this coordinate system use `sf::st_transform(county_one, crs = 6345)`. Afterwards you convert to `Spatial` and then `owin`. Note choosing the relevant projection is a science, and I don't know much about it, so do a bit of research if you want to make sure you don't get too distorted results. — Ege Rubak, Jan 06 '20 at 06:38

score 1 · Answer 2 · answered Oct 01 '20 at 15:22

just want to note here that coercion methods for sf classes are now registered (if that's the right word) by the sf package. I don't fully understand the R magic that finds methods, but it does work.

library(sf)
library(spatstat)

> methods(as.owin)
 [1] as.owin.boxx                    as.owin.data.frame              as.owin.default                 as.owin.distfun                
 [5] as.owin.dppm                    as.owin.funxy                   as.owin.im                      as.owin.influence.ppm          
 [9] as.owin.kppm                    as.owin.layered                 as.owin.leverage.ppm            as.owin.linfun                 
[13] as.owin.linnet                  as.owin.lintess                 as.owin.lpp                     as.owin.lppm                   
[17] as.owin.msr                     as.owin.MULTIPOLYGON*           as.owin.nnfun                   as.owin.owin                   
[21] as.owin.POLYGON*                as.owin.ppm                     as.owin.ppp                     as.owin.psp                    
[25] as.owin.quad                    as.owin.quadratcount            as.owin.quadrattest             as.owin.rmhmodel               
[29] as.owin.sf*                     as.owin.sfc*                    as.owin.sfc_MULTIPOLYGON*       as.owin.sfc_POLYGON*           
[33] as.owin.SpatialGridDataFrame*   as.owin.SpatialPixelsDataFrame* as.owin.SpatialPolygons*        as.owin.tess                   
see '?methods' for accessing help and source code

So, assuming you have properly projected your data (as noted by @Ege Rubak), calling as.owin directly should work:

library(tigris)

county <- county_subdivisions(state = "Alabama", 
                              county = "Lee", 
                              class = "sf", 
                              progress_bar = FALSE)

county <- st_union(county)

county <- st_transform(county, crs = 6345)

window <- as.owin(county)

plot(window)

Use sf polygon object as window in spatstat

2 Answers2

Linked