[R-sig-Geo] Prevent minor discrepancy in pixel resolution using Spatstat

[Roger Bivand]

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