[R] Spatstat: xy binary data into mask type to use in owin(mask=)

[Rolf Turner]

In the first instance, questions about contributed packages should be
addressed to the package maintainers rather than to the R-help list.

On 20/10/2009, at 2:06 AM, Javier PB wrote:

> Dear users,
>
> I am trying to export polygons from Arcmap into Spatstat to run some
> simulations using functions available in Spatstat package.
>
> One particular area to be exported is formed by a number of polygons
> defining the external boundaries of the area (as a groups of  
> islands) and a
> number of polygons inside the previous ones,  as “holes” not to be
> considered as part of the area.
>
> I have managed to export one of these areas using owin(poly=list()),
> including in the list the number of external boundaries and holes,  
> but the
> number of polygons that constitute each area is large and I was  
> wondering if
> this could be easily done using owin(mask=my_area), using only one  
> single
> matrix file "my_area" with TRUE, FALSE values defining the whole area.
>
> Here is where I got stuck. The file I have created in ArcMap  
> contains three
> columns: x and y coordinates, and the values at points x[i], y[j]  
> (1 if the
> area is TRUE and 0 is the area is FALSE).
>
> How can I covert this xy file into a mask type "my_area" that I can  
> read
> using owin(mask=my_area)?

It's not clear to me that you can --- it depends on the nature of the  
points
specified by x[i] and y[j].  If these are equispaced then I think you  
can
do it.

Suppose that you have a data frame ``X'' with columns ``x'', ``y'', and
``value'', where value is logical.  Suppose also that the unique values
of x are equispaced on [a,b] and the unique values of y are equispaced
on [c,d].  Finally suppose that *all* possible pairs of (x[i],y[j])
coordinates appear in your data frame X, i.e. nrow(X) = m*n where
m = the number of unique values of x and n = the number of unique values
of y.  Then you can create a mask from X as follows:

# Begin
	m <- as.mask(xy=X[,c("x","y")])
	ij <- with(X,nearest.raster.point(x[!value],y[!value],m))
	m$m[do.call(cbind,ij)] <- FALSE
# End

(Note that even if X$value is left as a numeric vector of 1's and 0's,
``!value'' will coerce it to logical appropriately.)

This will give you a window of type mask with bounding box [a,b] x [c,d]
which will ``be TRUE'' at those pixel centers (x[i],y[j]) where the
corresponding value of ``value'' is TRUE, and FALSE otherwise.

E.g.

	x0 <- seq(1,2,length=40)
	y0 <- seq(2,6,length=50)
	X  <- expand.grid(x=x0,y=y0)
	set.seed(42)
	X$value <- sample(c(TRUE,FALSE),2000,TRUE,prob=c(0.8,0.2))

Then do the three lines between # Begin and # End.  Then do

	plot(m)

to make sure it looks as it should.

There may be a sexier way of accomplishing the task.  If there is,
perhaps Adrian will chip in.

HTH.

	cheers,

		Rolf Turner
######################################################################
Attention: 
This e-mail message is privileged and confidential. If you are not the 
intended recipient please delete the message and notify the sender. 
Any views or opinions presented are solely those of the author.

This e-mail has been scanned and cleared by MailMarshal 
www.marshalsoftware.com
######################################################################