[R-sig-Geo] adding mask to spline raster

[Pascal Title]

Hi,

I am creating interpolations of global marine diversity via thin plate
splines, using the interpolate() function in the raster package. The
problem I am having is that if I plot the spline interpolation and then
overlay the plot with a continents polygon, it becomes clear that certain
areas are interpolated as having high diversity where they shouldn't due to
the configurations of continents, and the interpolation continuing "under"
the continents.

Is there any way to provide a mask either as a raster or as a spatial
polygon to the interpolate function? If not, can someone suggest a
different R package or method that may be better suited to this situation?

Here, I provide a dummy example. Areas with data on the outside of a ring
lead to the inside of the ring having high interpolated values. But if we
imagine that this was marine diversity, and an inland lake in the center
with no data, then we would want the "ring continent" to influence the
interpolation and "mask" the areas within the continent.


require(raster)

require(fields)

require(rgeos)


#making a box

b1 <- readWKT("POLYGON((200 200,300 200, 300 100, 200 100, 200 200))")

b2 <- readWKT("POLYGON((100 300, 400 300, 400 0, 100 0, 100 300))")

box <- gDifference(b2, b1)


#generate samples

b3 <- gDifference(gBuffer(b2, width=50), b2)

pts <- spsample(b3, 50, type='random')


b4 <- gDifference(gBuffer(b3, width=50), gEnvelope(pts))

pts2 <- spsample(b4, 50, type='random')


#create raster and populate

r <-raster(xmn=-100,xmx=600,ymn=-200,ymx=500,res=c(10,10))

r <- rasterize(gUnion(pts,pts2), r)

cells <- cellFromPolygon(r, b3)[[1]]

cells <- cells[which(!is.na(r[cells]))]

r[cells] <- round(rnorm(length(cells), mean=500, sd=20))


cells <- cellFromPolygon(r, b4)[[1]]

cells <- cells[which(!is.na(r[cells]))]

r[cells] <- round(rnorm(length(cells), mean=50, sd=20))


#generate thin plate spline

xy <- data.frame(xyFromCell(r, 1:ncell(r)))

tps <- Tps(xy, values(r))

fine <- disaggregate(r, 2)

p <- raster(fine)

p <- interpolate(p, tps)


#plot

colors <- colorRampPalette(c('blue','yellow','red'))

plot(p, col=colors(100))

plot(box,add=TRUE, col='black')

plot(r, add=TRUE)



Thanks!
-Pascal


-- 
Pascal Title
PhD candidate | Rabosky Lab
<http://www-personal.umich.edu/~drabosky/Home.html>
Dept of Ecology and Evolutionary Biology
University of Michigan
ptitle at umich.edu

pascaltitle.weebly.com

	[[alternative HTML version deleted]]