spatstat <--> sf conversions

[edzer]

I tried to implement a number of sf -> spatstat conversions, and vice versa. A review by the spatstat team e.g. @rubak would be appreciated. I'm in particular not sure whether the observation windows for the psp objects make sense.

Here is a demo:

suppressPackageStartupMessages(library(spatstat))
suppressPackageStartupMessages(library(sf))

# png("/tmp/spa%03d.png")

p1 = st_point(0:1)
p2 = st_point(1:2)
p3 = st_point(c(-1,2))
p = st_sfc(p1, p2, p3)
as.ppp(p)
# Planar point pattern: 3 points
# window: rectangle = [-1, 1] x [1, 2] units
try(as.ppp(st_set_crs(p, 4326)))
# Error in as.ppp.sfc(st_set_crs(p, 4326)) : 
#   Only projected coordinates may be converted to spatstat class objects

sf = st_sf(geom = p)
try(as.ppp(sf))
# Planar point pattern: 3 points
# window: rectangle = [-1, 1] x [1, 2] units
sf = st_sf(a = 1:3, geom = p)
as.ppp(sf)
# Marked planar point pattern: 3 points
# marks are numeric, of storage type  ‘integer’
# window: rectangle = [-1, 1] x [1, 2] units
sf = st_sf(a = 1:3, b=3:1, geom = p)
as.ppp(sf) # warns
# Marked planar point pattern: 3 points
# marks are numeric, of storage type  ‘integer’
# window: rectangle = [-1, 1] x [1, 2] units
# Warning message:
# In as.ppp.sf(sf) : only first attribute column is used for marks

w = st_as_sfc(st_bbox(st_sfc(p1, p2)))
sf = st_sf(a = 1:3, geom = p)
(p0 = rbind(st_sf(a = 0, geom = w), sf))
# Simple feature collection with 4 features and 1 field
# geometry type:  GEOMETRY
# dimension:      XY
# bbox:           xmin: -1 ymin: 1 xmax: 1 ymax: 2
# epsg (SRID):    NA
# proj4string:    NA
#   a                           geom
# 1 0 POLYGON ((0 1, 1 1, 1 2, 0 ...
# 2 1                    POINT (0 1)
# 3 2                    POINT (1 2)
# 4 3                   POINT (-1 2)
try(as.ppp(p0)) # errors: one point outside window
# Error in `marks<-.ppp`(`*tmp*`, value = value) : 
#   number of rows of data frame != number of points
# In addition: Warning message:
# 1 point was rejected as lying outside the specified window 

w = st_as_sfc(st_bbox(p))
sf = st_sf(a = 1:3, geom = p)
(p0 = rbind(st_sf(a = 0, geom = w), sf))
# Simple feature collection with 4 features and 1 field
# geometry type:  GEOMETRY
# dimension:      XY
# bbox:           xmin: -1 ymin: 1 xmax: 1 ymax: 2
# epsg (SRID):    NA
# proj4string:    NA
#   a                           geom
# 1 0 POLYGON ((-1 1, 1 1, 1 2, -...
# 2 1                    POINT (0 1)
# 3 2                    POINT (1 2)
# 4 3                   POINT (-1 2)
as.ppp(p0)
# Marked planar point pattern: 3 points
# marks are numeric, of storage type  ‘double’
# window: polygonal boundary
# enclosing rectangle: [-1, 1] x [1, 2] units

library(stars)
# Loading required package: abind
tif = system.file("tif/L7_ETMs.tif", package = "stars")
s = adrop(read_stars(tif)[,,,1]) > 70
plot(s)

m = as.owin(s)
plot(m)

table(m$m)

# FALSE  TRUE 
# 39494 83354 

# as.owin.sf, as.owin.sfc_*
nc = st_read(system.file("gpkg/nc.gpkg", package="sf"))
# Reading layer `nc.gpkg' from data source `/home/edzer/R/x86_64-pc-linux-gnu-library/3.6/sf/gpkg/nc.gpkg' using driver `GPKG'
# Simple feature collection with 100 features and 14 fields
# geometry type:  MULTIPOLYGON
# dimension:      XY
# bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
# epsg (SRID):    4267
# proj4string:    +proj=longlat +datum=NAD27 +no_defs
try(as.owin(nc)) # should be projected
# Error in as.owin.sfc_MULTIPOLYGON(st_geometry(W), ...) : 
#   Only projected coordinates may be converted to spatstat class objects
nc = st_transform(nc, 32119)
plot(as.owin(nc), col = 'grey')

plot(as.owin(st_geometry(nc)), col = 'grey')

sq = rbind(c(-1,-1), c(1, -1), c(1,1), c(-1,1), c(-1,-1))
pol = st_polygon(list(0.5 * sq, sq[5:1,] * 0.45)) # w hole
plot(as.owin(pol), col = 'grey')

plot(as.owin(st_sfc(pol)), col = 'grey')

mpol = st_multipolygon(list(
	list(sq, sq[5:1,] * 0.9),
	list(sq * 2, sq[5:1,] * 1.8)))
plot(as.owin(mpol), col = 'grey')

plot(as.owin(st_sfc(mpol)), col = 'grey')

plot(as.owin(st_sfc(pol, mpol)), col = 'grey')

plot(as.owin(st_sf(a=1:2, st_sfc(pol, mpol))), col = 'grey')

o = as.owin(st_sf(a=1:2, st_sfc(pol, mpol)))

st_as_sfc(o)[[1]]
# POLYGON ((2 2, -2 2, -2 -2, 2 -2, 2 2), (-1.8 -1.8, -1.8 1.8, 1.8 1.8, 1.8 -1.8, -1.8 -1.8))
st_as_sfc(o)[[2]]
# POLYGON ((1 1, -1 1, -1 -1, 1 -1, 1 1), (-0.9 -0.9, -0.9 0.9, 0.9 0.9, 0.9 -0.9, -0.9 -0.9))
st_as_sfc(o)[[3]]
# POLYGON ((0.5 0.5, -0.5 0.5, -0.5 -0.5, 0.5 -0.5, 0.5 0.5), (-0.45 -0.45, -0.45 0.45, 0.45 0.45, 0.45 -0.45, -0.45 -0.45))

plot(st_as_sfc(o), col = 'blue', main = 'st_as_sfc(o)')

plot(st_as_sf(o), col = 'blue', main = 'st_as_sf(o)')

data(nztrees)
qNZ <- quadratcount(nztrees, nx=4, ny=3)
ts = as.tess(qNZ)
plot(st_as_sfc(ts))

ls = st_linestring(rbind(c(0,0), c(1,1), c(2,0)))
plot(as.psp(ls))

mls = st_multilinestring(list(rbind(c(0,0), c(1,1), c(2,0)), rbind(c(3,3), c(4,2))))
plot(as.psp(mls))

plot(as.psp(st_sfc(ls)))

plot(as.psp(st_sfc(mls)))

plot(as.psp(st_sfc(ls, mls)))

[edzer]

Note that this is done in a branch called "ppp".

[rubak]

This looks great @edzer! I will try to have a closer look this coming week.

[gueyenono]

EDIT:

I installed the development version of the {stars} package and it worked

end EDIT.

I am not able to replicate this code chunk:

suppressPackageStartupMessages(library(spatstat))
suppressPackageStartupMessages(library(sf)) # r-spatial/sf@ppp (ppp branch of the github repo)
library(stars)

tif = system.file("tif/L7_ETMs.tif", package = "stars")
s = adrop(read_stars(tif)[,,,1]) > 70
plot(s)
m = as.owin(s) # does not work
plot(m) # does not work (because previous line of code)

I get this error when I run m = as.owin(s): Error in as.owin(s) : could not find function "as.owin" . I installed the "ppp" branch of {sf} and the development version of {spatstat}.

[gueyenono]

I installed the ppp branch of the {sf} package; however, I am not able to run this code: st_read(dsn = system.file("gpkg/nc.gpkg", package="sf")). I get the following error:

Error in .Call("_sf_CPL_read_ogr", PACKAGE = "sf", datasource, layer,  : 
  Incorrect number of arguments (13), expecting 12 for '_sf_CPL_read_ogr'

[edzer]

@gueyenono for the stars problem, you need to install the github version of stars, branch master; for the nc reading problem: maybe remove sf altogether (by hand) and then reinstall it. If it persists, please give your devtools::sessions_info() output.

[gueyenono]

@gueyenono for the stars problem, you need to install the github version of stars, branch master; for the nc reading problem: maybe remove sf altogether (by hand) and then reinstall it. If it persists, please give your devtools::sessions_info() output.

Removing sf by hand and reinstalling it worked. Thank you.

[edzer]

I merged this branch now into master. Have you had time to look at it, @rubak ?

[rubak]

I have only started looking, but I have an example where a ppp with a polygonal owin with holes is mis-intepreted. Funny thing is that direct conversion of the owin works fine. Reproducible example:

suppressPackageStartupMessages(library(spatstat))
suppressPackageStartupMessages(library(sf))

W <- Window(gordon)
plot(W)

A point in a polygonal hole is correctly classified as “outside” by spatstat:

inside.owin(4, -27, W)
#> [1] FALSE

Converting the owin directly to sf makes it into a single polygon with holes and
things work correctly:

W_sf <- st_as_sf(W)
W_sf
#> Simple feature collection with 1 feature and 0 fields
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: -26.40848 ymin: -36.32095 xmax: 26.40848 ymax: 36.32095
#> epsg (SRID):    NA
#> proj4string:    NA
#>                             geom
#> 1 POLYGON ((22.0448 -14.20711...
p_sf <- st_point(matrix(c(4,-27),ncol=2))
st_contains(W_sf, p_sf) # Point not contained in polygon (correct)
#> Sparse geometry binary predicate list of length 1, where the predicate was `contains'
#>  1: (empty)

However if W is the window of a ppp (in this case empty, but same issue with points)
it is converted incorrectly:

X <- ppp(numeric(0), numeric(0), window = W)
W_sf2 <- st_as_sf(X)
W_sf2
#> Simple feature collection with 1 feature and 1 field
#> geometry type:  MULTIPOLYGON
#> dimension:      XY
#> bbox:           xmin: -26.40848 ymin: -36.32095 xmax: 26.40848 ymax: 36.32095
#> epsg (SRID):    NA
#> proj4string:    NA
#>    label                           geom
#> 1 window MULTIPOLYGON (((22.49851 -1...
st_contains(W_sf2, p_sf) # Point allegedly in polygon (wrong)
#> Sparse geometry binary predicate list of length 1, where the predicate was `contains'
#>  1: 1

Probably, the problem is due to the conversion of the polygon through edges here:

sf/R/spatstat.R

Lines 1 to 9 in 1e095d4

	window_polygons_from_edges = function (w) {
	mw = as.matrix(w$ends)
	lst1 = lapply(seq_len(NROW(mw)), function(i) st_linestring(matrix(mw[i,], 2, byrow = TRUE)))
	p0 = st_polygonize(do.call(c, do.call(st_sfc, lst1)))
	if (length(p0) > 1) # multiple POLYGONs, returned as sfc_
	do.call(c, st_collection_extract(p0, "POLYGON")) # MULTIPOLYGON
	else
	st_cast(p0, "POLYGON")
	}

I assume the same problem is present when converting psp to sf. Maybe the conversion of owin to polygon could be inspired by maptools:
https://github.com/cran/maptools/blob/14a028b493b22969d03af5ed221b5fe992cc1b67/R/spatstat1.R#L69-L97

[rubak]

Just checked that the problem is the same for psp.

Now that I looked further into the code I see that the maptools code for polygons is already in st_as_sfc.owin, so we just need to use that code in st_as_sf.ppp and st_as_sf.psp, so the fix should be simple.

I also noticed there are some issues with multiple columns of marks (at least for psp -- seems OK for ppp)...

[edzer]

Thanks! I think st_as_sf.ppp and st_as_sf.psp now set the right window.

[bastistician]

Should check_ring_dir (or similar) also be called in as.owin.POLYGON and as.owin.MULTIPOLYGON? Currently the following fails:

nc <- st_read(system.file("shape/nc.shp", package="sf"), quiet = TRUE)
spatstat::as.owin(nc$geometry[[1]])
#> Error in spatstat::owin(bb[c("xmin", "xmax")], bb[c("ymin", "ymax")],  : 
#>   Area of window is negative;
#>  check that all polygons were traversed in the right direction

[edzer]

Yes, that makes sense!