Prevent minor discrepancy in pixel resolution using Spatstat
On Wed, 17 Feb 2010, Tyler Dean Rudolph wrote:
When I specify eps=25 (m resolution) in a call to distmap...
library(spatstat)
temp <- data.frame(x=-670049.2, y=752814.6, X2=-669961.9, Y2=752648.2)
subwin<-owin(xrange=c(-100, 100) + sort(c(temp$x, temp$X2)),
yrange=c(-100, 100) + sort(c(temp$y, temp$Y2)))
trajpsp<-psp(temp$x,temp$y, temp$X2, temp$Y2, subwin, check=TRUE)
plot(trajmap<-distmap(trajpsp, eps=25))
trajmap$xstep
[1] 23.94167
trajmap$ystep
[1] 24.42667 ...The resulting resolution is slightly different than 25x25 metres.
I think that you are seeing the impact of your input window not being
exactly divisible by your eps in as.mask() and then getting nr and nc that
do not partition the x and y ranges into discrete 25m blocks, because of
ceiling().
if (!is.null(eps)) {
nr <- ceiling(diff(w$yrange)/eps)
nc <- ceiling(diff(w$xrange)/eps)
}
Browse[2]> nr
[1] 15
Browse[2]> nc
[1] 12
With:
subwin<-owin(xrange=c(-100, 110+2.7) + sort(c(temp$x, temp$X2)),
yrange=c(-100, 90+18.6) + sort(c(temp$y, temp$Y2)))
you get something like what you want, with a lot of fiddling. Have you
considered not setting eps, using spatstat.options("npixel") to set nr and
nc, and tailoring the window to suit 25m (yrange == eps*nr, etc.). If you
know the xrange, yrange of the raster you want to compare with, maybe use
that as the source for the window, spatstat.options("npixel") to feed
through nr and nc, then you could fix eps indirectly?
Hope this helps,
Roger
So now when I use it as a mask:
trajmap$v<-ifelse(trajmap$v>50, 0, 1) library(maptools)
...and convert it to a RasterLayer for matrix algebra with another raster
trajrast<-raster(as.SpatialGridDataFrame.im(from=trajmap)) tempmap<-crop(habmap, extent(trajrast)) # tempmap.grd is attached top<-table(getValues(tempmap*trajrast))
Error in compare(c(e1, e2)) : Different bounding box Error in getValues(tempmap * trajrast) : error in evaluating the argument 'x' in selecting a method for function 'getValues' ...it doesn't work because the resolution is not exactly equal. Now I can sometimes trick it when the number of rows and cells are still equal, BUT... extent(trajrast)<-extent(tempmap)
top<-table(getValues(tempmap*trajrast))
Error in compare(c(e1, e2)) : nrows different Error in getValues(tempmap * trajrast) : error in evaluating the argument 'x' in selecting a method for function 'getValues' ...This is not always the case! How can I prevent the numbers of rows and columns from being so slightly different because of a minor different in resolution? Otherwise put, how can I make the actual resolution of my distmap image exactly 25m as desired? Tyler
Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: Roger.Bivand at nhh.no