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@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
p...@stat.auckland.ac.nz
http://www.stat.auckland.ac.nz/~paul/
______________________________________________
R-help@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.
