R help - Jul 2013 - Masking oceans using polypath

Hi R-help

I am trying to mask the ocean from an image plot I have made.


Here is some example code:

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


As you can see, the polypath is on both sides of the country outline. Any
help is much appreciated.


Louise

	[[alternative HTML version deleted]]

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")
-- 
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/