[R-sig-Geo] spatstat crashes R when converting tess to sp objects

[O'Hanlon, Simon J]

Dear list,
I am experiencing a reoccurring crash in R when trying to convert a tess object from spatstat into an sp object. It seems to occur only when I create a tess object using the dirichlet function, specifically when specifying a polygonal owin object in the window argument. I am using R 2.15.1 on Mac OS X 10.7.5. (Platform: x86_64-apple-darwin9.8.0/x86_64 (64-bit) ). I have spatstat ver. 1.29-0 and sp 1.0.1. This also happened on my Windows 7 64-bit box (albeit without the helpful error message), with sp 0.9.99, spatstat 1.28.1, and R version 2.15.1 (2012-06-22) Platform: x86_64-pc-mingw32/x64 (64-bit).

I am trying to create voronoi/thiessen polygons from a point pattern, limited to within country borders.

The error thrown up is this:

*** caught segfault ***
address 0x0, cause 'memory not mapped'


Traceback:
 1: .Call("Polygons_c", srl, ID, PACKAGE = "sp")
 2: Polygons(pieces, id)
 3: owin2Polygons(y[[i]], nam[i])
 4: tess2SP(v)

Here is a small reproducible example that you should be able to copy and paste. BUT PLEASE NOTE THIS CRASHES MY CURRENT R SESSION WHEN I TRY TO RUN IT. I'd hate for anyone to lose some work! Can anyone else recreate my crash? I am guessing that the cause of my problems lies within the Polygons_c file within spatstat?


#Reproducible example

#Required libraries
library(sp); library(spatstat); library(maptools); library(rgeos)


# Data - two columns with x and y coords
p <- as.data.frame( structure(c(36.8547, 8.66667, 38.9782, 34.3265, 10.0439, 15.3166, 
15.4833, 32.8404, -6.1833, 17.8167, -16.3522, 3.46667, -1.5489, 
-0.2818, 12.449, -4.1333, -17.7167, -16.8231, 14.066, -20.5833
), .Dim = c(10L, 2L), .Dimnames = list(NULL, c("longitude", "latitude"
))))

#	Download and unzip country borders data
oldwd <- getwd(); tmp <- tempdir(); setwd(tmp)
url <- "http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/50m/cultural/50m-admin-0-countries.zip"
dest <- paste(tmp,"\\tmp.zip",sep="")
download.file(url,dest)	#File is 1.3Mb
unzip(dest)

# Read in world shapefile
wld <- readShapePoly("ne_50m_admin_0_countries")

#Subset to the countries just in Africa
countries <- wld at data[order(wld at data$NAME),]
afr <- wld[ match( countries$NAME[c(4,7,25,30,36,37,39,41,43,44,48,49,50,59,64,66,67,69,77,78,82,89,90,107,113,121,122,123,130,133,136,137,145,146,151,158,159,178,181,186,188,190,191,197,198,199,207,209,215,217,220,225,236,240,241)] , wld at data$NAME ) , ]
#Melt into one polygon
africa <- gUnaryUnion( afr )

#Create polygonal observation window for voronoi polygons
win <- as.owin(africa)
points <- ppp( p[,1] , y = p[,2] , window = win )
# Create voronoi polygons
v <- dirichlet( points )
# Does not seem to work on my system
vp <- as( v , "SpatialPolygons" )
# End of example


Thanks in advance for any help,


Simon

------------------------------------
Simon O’Hanlon
Postgraduate Researcher
Department of Infectious Disease Epidemiology
Imperial College London