[R] Masking oceans using polypath

[Paul Murrell]

Hi

The outline for France from that database is a little bit pesky because 
it consists of multiple segments that travel in different directions 
around the French border.  Here is some (possibly jetlagged) code that 
puts the segments all in the same direction and produces a nicer result ...

breaks <- which(is.na(outline$x))
starts <- c(1, breaks + 1)
ends <- c(breaks - 1, length(outline$x))
sections <- mapply(":", starts, ends)
xx <- lapply(sections, function(i) { outline$x[i] })
yy <- lapply(sections, function(i) { outline$y[i] })
for (i in 2:length(xx)) {
     if (xx[[i - 1]][length(xx[[i - 1]])] != xx[[i]][1] ||
         yy[[i - 1]][length(xx[[i - 1]])] != yy[[i]][1]) {
         xx[[i]] <- rev(xx[[i]])
         yy[[i]] <- rev(yy[[i]])
     }
}
image(x=seq(from=xbox[1], to=xbox[2], length=10),
        y =seq(from=ybox[1], to=ybox[2], length=10),
        z = outer(1:10, 1:10, "+"),
        xlab = "lon", ylab = "lat", xlim=xbox, ylim=ybox)
polypath(c(unlist(xx), NA, c(rev(xbox), xbox)),
           c(unlist(yy), NA, rep(ybox, each=2)),
           col="light blue", rule="evenodd")

... Hope that works for you too.

Paul

On 08/06/13 04:02, Marc Girondot wrote:
> I like very much the solution proposed by Paul Murrell to mask the ocean
> (in fact, all that is not the country displayed) in a plot. It works
> pretty well for Australia. However, when I try to apply the same code
> for France, it fails. I search for the reason but I can't find. Here is
> the code for Australia, after the same code for France with the display
> problem.
> Thanks a lot for any advice on the reason for the discrepancy.
>
> Marc Girondot
>
>
> ### For Australia
> library(mapdata)
>
> outline <- map("worldHires", regions="Australia", exact=TRUE,
>                  plot=FALSE) # returns a list of x/y coords
>
> xrange <- range(outline$x, na.rm=TRUE) # get bounding box
> yrange <- range(outline$y, na.rm=TRUE)
>
> xbox <- round(xrange + c(-2, 2))
> ybox <- round(yrange + c(-2, 2))
>
> image(x=seq(from=xbox[1], to=xbox[2], length=10),
>         y =seq(from=ybox[1], to=ybox[2], length=10),
>         z = outer(1:10, 1:10, "+"),
>         xlab = "lon", ylab = "lat", xlim=xbox, ylim=ybox)
>
>
> # create the grid path in the current device
>
> subset <- !is.na(outline$x)
> polypath(c(outline$x[subset], NA, c(rev(xbox), xbox)),
>            c(outline$y[subset], NA, rep(ybox, each=2)),
>            col="light blue", rule="evenodd")
>
> ## Now here is the same code for France... not so beautiful.
>
> ### For France
> library(mapdata)
>
> outline <- map("worldHires", regions="France", exact=TRUE,
>                  plot=FALSE) # returns a list of x/y coords
>
> xrange <- range(outline$x, na.rm=TRUE) # get bounding box
> yrange <- range(outline$y, na.rm=TRUE)
>
> xbox <- round(xrange + c(-2, 2))
> ybox <- round(yrange + c(-2, 2))
>
> image(x=seq(from=xbox[1], to=xbox[2], length=10),
>         y =seq(from=ybox[1], to=ybox[2], length=10),
>         z = outer(1:10, 1:10, "+"),
>         xlab = "lon", ylab = "lat", xlim=xbox, ylim=ybox)
>
>
> # create the grid path in the current device
>
> subset <- !is.na(outline$x)
> polypath(c(outline$x[subset], NA, c(rev(xbox), xbox)),
>            c(outline$y[subset], NA, rep(ybox, each=2)),
>            col="light blue", rule="evenodd")
>
>
>
>
> Le 16/07/13 23:42, Paul Murrell a écrit :
>> Hi
>>
>> There are a couple of problems:
>>
>> 1.
>> Your 'outline' is much bigger than it needs to be.  For example, the
>> following produces just Australia ...
>>
>> outline <- map("worldHires", regions="Australia", exact=TRUE,
>>                 plot=FALSE) # returns a list of x/y coords
>>
>> 2.
>> The outline you get still has NA values in it (the coastline of
>> Australia in this database is a series of distinct lines;  I don't
>> know why that is).  Fortunately, you can stitch Australia back
>> together again like this ...
>>
>> subset <- !is.na(outline$x)
>> polypath(c(outline$x[subset], NA, c(xbox, rev(xbox))),
>>           c(outline$y[subset], NA, rep(ybox, each=2)),
>>           col="light blue", rule="evenodd")
>>
>> With those two changes, I get a much better result.  You may want to
>> fiddle a bit more to add Tasmania and bits of PNG and Indonesia back
>> in the picture OR you could replace these problems with a different
>> approach and use a more sophisticated map outline via a "shapefile"
>> (see the 'sp' package and the 'maptools' package and possibly the
>> R-sig-geo mailing list).
>>
>> Hope that helps.
>>
>> Paul
>>
>> On 07/16/13 23:17, Louise Wilson wrote:
>>> library(mapdata)
>>>
>>> image(x=110:155, y =-40:-10, z = outer(1:45, 1:30, "+"),
>>>
>>>         xlab = "lon", ylab = "lat")
>>>
>>> outline <- map("worldHires", plot=FALSE) # returns a list of x/y coords
>>>
>>> xrange <- range(outline$x, na.rm=TRUE) # get bounding box
>>>
>>> yrange <- range(outline$y, na.rm=TRUE)
>>>
>>> xbox <- xrange + c(-2, 2)
>>>
>>> ybox <- yrange + c(-2, 2)
>>>
>>> # create the grid path in the current device
>>>
>>> polypath(c(outline$x, NA, c(xbox, rev(xbox))),
>>>
>>>        c(outline$y, NA, rep(ybox, each=2)),
>>>
>>>        col="light blue", rule="evenodd")
>>
>
>
>
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>

-- 
Dr Paul Murrell
Department of Statistics
The University of Auckland
Private Bag 92019
Auckland
New Zealand
64 9 3737599 x85392
paul at stat.auckland.ac.nz
http://www.stat.auckland.ac.nz/~paul/