[R-sig-Geo] How to mask out

[Roger Bivand]

On Tue, 4 Jan 2005, Edzer J. Pebesma wrote:

> José, you would have to do the following:
> 1. obtain the polygons the make up the Iberian peninsula
> 2. overlay the interpolation grid with the polygons:
>    for each point on the grid, determine whether it is
>    inside or outside the peninsula;
> 3. put these (logical) values in a matrix of the same
>    form of the z argument in the list returned by interp.new
>    (the matrix in the form that "image" likes)
> 4. using this logical matrix, say B, mask out the "outside"
>    values in the interpolated matrix, say A$z, by
>    > A$z[B] = NA
> 
> I don't know whether maps can give you the complete polygons,
> I think they work from topology; if rings are formed on the
> fly, I don't know how to get them. Try:
> 
> library(maps)
> x = map("world", "spain")
> polygon(cbind(x$x,x$y),col="red")
> 
> and you'll see that x does not contain one single polygon
> for the big peninsula, but many fragments.
> 
> Points in polygon are not so hard to get in R.
> My guess is that Roger might be able to help you further.

Well, after that introduction, I'm going to have to!

If continental Spain is enough, then:

library(maps)
library(mapdata)
x <- map("worldHires", "Spain", exact=TRUE, fill=TRUE, plot=FALSE)

will extract just the polygon (with the arcs in the correct direction - 
fill=TRUE is needed for that, otherwise you get line segments in 
arbitrary directions). Then:

library(splancs)
grid <- gridpts(as.points(x), xs=0.1, ys=0.1)
plot.new()
plot.window(xlim=range(x$x), ylim=range(x$y))
polygon(x)
points(grid, pch=".")

gives a grid at 0.1 degree spacings. My guess would be that you actually 
need a projected polygon, and that you can generate your grid from your 
interp() output x and y components, and use the splancs package inout() 
to give a logical vector, having set extrap=TRUE in interp(). You will 
need to be careful overlaying the mask from inout() on the interpolated 
values, but it can be done. Unfortunately, mapproject() in package mapproj 
returns coordinates suited to plotting but not other use, so you may need 
to write the lat-lon coordinates of the polygon out, and use proj to 
project before reading them back in again.

Hope this helps,

Roger



> --
> Edzer
> 
> José Agustín García García wrote:
> > Hello
> > 
> > I am traying to interpolate several rainfalls data over the Iberian 
> > Peninsula and display the results over a map of the Peninsula. I have 
> > used the akima package to make the interpolation and the maps package to 
> > draw the map. My cuestion is whow may I  mask out the output of the 
> > interp.new function (with extrapolation) with the iberian Peninsula map 
> > provided by the map package in order to restrict the image to the 
> > Iberian Peninsula ? . Any suggestions are wellcome
> > 
> > Thank in advance for your help
> > 
> > Agustin
> 
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Breiviksveien 40, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 93 93
e-mail: Roger.Bivand at nhh.no