Re: [R] Masking oceans using polypath
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
Re: [R] Masking oceans using polypath
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")
>
--
__
Marc Girondot, Pr
Laboratoire Ecologie, Systématique et Evolution
Equipe de Conservation des Populations et des Communautés
CNRS, AgroParisTech et Université Paris-Sud 11 , UMR 8079
Bâtiment 362
91405 Orsay Cedex, France
Tel: 33 1 (0)1.69.15.72.30 Fax: 33 1 (0)1.69.15.73.53
e-mail: marc.giron...@u-psud.fr
Web: http://www.ese.u-psud.fr/epc/conservation/Marc.html
Skype: girondot
[[alternative HTML version deleted]]
__
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.
Re: [R] Masking oceans using polypath
Many thanks, this works beautifully.
Louise
On Wed, Jul 17, 2013 at 7:42 AM, Paul Murrell wrote:
> 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
> p...@stat.auckland.ac.nz
> http://www.stat.auckland.ac.**nz/~paul/
Re: [R] Masking oceans using polypath
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
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.
[R] 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]]
__
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.
