Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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.

  require(spatstat)
  data(lansing)
  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.

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

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

share|improve this question
    
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
add comment

1 Answer

up vote 1 down vote accepted

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)
}
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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