[R-sig-Geo] overlay fx, consistent topology, sample a point on a 'patch edge'

[Roger Bivand]

On Thu, 7 Feb 2008, Edzer Pebesma wrote:

> Ilona Naujokaitis-Lewis wrote:
>> NO
>>
>> Dear list-serve,
>> Thanks in advance to all those who help out with the inquiries, it is has
>> helped me numerous times. Here go my questions...
>>
>> I am trying to vary existing landscapes which are composed of habitat
>> (patches) and non-habitat. My goals are to vary the number of patches in a
>> landscape, and the size of each patch.
>>
>> My landscape files are originally raster files, which I have converted to
>> ascii format. When importing into R, they are of class SpatialGridDataFrame,
>> are fully gridded. Each cell is represented by a 0: non-habitat, or a value
>> ranging between >0 and <=1, which represents varying degrees of habitat
>> quality. Patches are simply identified where adjoining cells are
>>
>>> 0.
>>>
>>
>> The approach I have taken is to create a mask of my ascii landscape file so
>> that already existing patches are masked. This allows me to sample a point
>> in the region of my landscape that is identified as non-habitat. I then
>> apply the dilate function to essentially 'grow' a patch.
>>
>> I have followed the code from the following link:
>> http://finzi.psych.upenn.edu/R/Rhelp02a/archive/115868.html
>> The only difference is that once I have effected the dilate function, which
>> is of class 'owin', I then:
>> 1) coerce "owin" to "im" object class
>> 2) coerce "im" object class to a SpatialGridDataframe
>>
>> My problems are:
>> 1. How to overlay my final 2 layers:  1) dilated object (object class
>> SpatialGridDataFrame) and 2) original landscape patch layer (also object
>> class SpatialGridDataFrame). I need the resulting object to be an object of
>> class SpatialGridDataFrame with only 1 band consisting of original patches,
>> the new patches, and non-habitat so I can perform additional data
>> manipulations. The overlay function does not appear to work with 2 objects
>> of this class. Any alternative suggestions?
>>
> Method overlay only combines objects of different class. If you have two
> grids and want to create a third, and all three have the same
> topology/number of cells, then simply use vectorized expressions like
>
> grd$out = grd$dilated != grd$original

This is possible, but dilute.owin() not only extends the im object where 
the dilated region would cross the existing boundary, but also clips off 
parts not containing dilated cells. In this modified version:

dilate1.owin <- function (w, r, ..., include.bb=FALSE)
{
     verifyclass(w, "owin")
     if (r < 0)
         stop("r must be nonnegative")
     if (r == 0)
         return(w)
     if (w$type != "mask")
         w <- as.mask(w, ...)
     epsilon <- sqrt(w$xstep^2 + w$ystep^2)
     r <- max(r, epsilon)
     if (!include.bb) bb <- bounding.box(w)
     else bb <- w
     newbox <- owin(bb$xrange + c(-r, r), bb$yrange + c(-r, r))
     w <- rebound.owin(w, newbox)
     distant <- distmap(w)
     dil <- w
     dil$m <- (distant$v <= r)
     return(dil)
}

the input ranges are used, buffered out by the dilation range in each 
direction. So from the cited example:

> source("dilate1.owin")
> pts_owin_0.03 <- dilate1.owin(pts_owin, r=0.03, include.bb=TRUE)
> pts_owin_0.03
window: binary image mask
106 x 106 pixel array (ny, nx)
enclosing rectangle: [-0.03, 1.03] x [-0.03, 1.03] units
> x <- as.im(pts_owin_0.03)
> xx <- as(x, "SpatialGridDataFrame")
> xxx <- overlay(as(SGDF, "SpatialPixelsDataFrame"), xx)
> xx$orig[!is.na(xxx)] <- TRUE
> image(xx, "v")
> image(xx, "orig", col="black", add=TRUE)

gets you there more or less.

Roger

PS. The other questions are still open, the edge detection feels a bit 
like image analysis, say the biOps package for a range of filters to 
choose the cells on the edge of the patch (perhaps, untried).

>> 2. My dilated points object (SGDF) does not have topology identical to the
>> original layer - I think I need them to be the same for the overlay process,
>> and if yes, how do I change this. The cell size remains the same but the
>> offset and cell dimension are not identical.
>> For example:
>> > getGridTopology(dilate_SGDF)
>>   	                  X1  X2
>> 	cellcentre.offset 1050 890
>> 	cellsize            20  20
>> 	cells.dim           15  11
>> > getGridTopology(original_layer)
>>         	             x   y
>> 	cellcentre.offset 1030 870
>> 	cellsize            20  20
>> 	cells.dim           16  12
>>
> You could coerce the "wrong" one into a SpatialPointsDataFrame and the
> use function (method) overlay with the good grid to get the right indexes.
> --
> Edzer
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
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