i am using the very good spatial statistic package spatstat to produce probability maps of spatial point patterns. Here is a simple working result from the examples to try out.

  z <- relrisk(lansing)
  plot(z, main="Lansing Woods")

My question is rather simple (or maybe not). I produced a similar map as shown above and i want to crop it in R with a boundary shapefile. The spatstat package uses its own system (of windows and pixel images) so i am rather clueless. I know that i can format my shapefile to a so called mask with this command as.mask(as.owin(shapefile)), but i didn't get beyond this point.
Anyone knows how i can crop the produced image with this mask?

EDIT: Ahh, sry. I already found a solution around it, which works for me. I simply have to make a subset out of the spatial points before converting it to a ppp object.

st <- readOGR("GIS-Data","boundary_shape") 
wp <- spp[-which(is.na(over(spp,st)[,1])),] # spp is my spatial point shape
#Then use this command from the maptools package
plp <- as.ppp.SpatialPointsDataFrame(wp)

Nevertheless it would still be nice to know if it is possible to crop a spatstat pixel image.

  • I'm a bit confused by what you really want - do you want to know how to apriori specify a window for the point pattern, or do you want to know how to simply crop the plot to a pre-specified x-y domain? – Andy W Mar 3 '13 at 19:41
  • nja, i produced an image using the relrisk function and i wanted to show only a subset of the image to some folks. Therefore i want to know how to simply crop the plot to a x-y domain (a boundary polygons extent). However removing some of the points outside the polygon didn't really changed the image looks, therefore i will go with that for now. – Curlew Mar 3 '13 at 22:21

The best approach is to coerce your spatstat "image" object to an sp object and use the over function or clip in the raster package. Here is a function to perform the coercion to an sp SpatialGridDataFrame.

as.SpatialGridDataFrame.im <- function(from){
    offset=c(from$xrange[1] + 0.5 * from$xstep, from$yrange[1] + 0.5 * from$ystep)
      cellsize = c(diff(from$xrange)/from$dim[2], diff(from$yrange)/from$dim[1])
dim = c(from$dim[2], from$dim[1])
  gt = GridTopology(offset, cellsize, dim)
    m = t(from$v[nrow(from$v):1, ])
  data = data.frame(v = as.vector(m))
SpatialGridDataFrame(gt, data)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.