[R-sig-Geo] plotting a polygon with holes

Bjarke Christensen Bjarke.Christensen at sydbank.dk
Thu Jun 11 15:02:28 CEST 2009


Hi,

I have a shapefile/SpatialPolygonsDataFrame of coastlines. Each polygon
corresponds to an island, and there are no holes. I can plot this so that
the islands are shaded by using
> plot(islands, col="gray")

What I now want, is to plot the same information so that the ocean is blue
and the islands are transparent. Something like:
>image(myDEM)
>plot(ocean, col="lightblue", add=T)

which I would hope would allow the DEM to be visible on the islands, but
not in the ocean.

The approach I have been taking does not seem to work, so I would like to
as for your suggestions. My approach is to make a new map, with one, big
rectangular polygon containing all the islands as holes.
I do this by making another shapefile, containing one polygon which
consists of the (slightly expanded) bounding box of the first shapefile. I
then import both shapefiles into GRASS, run

       v.overlay ainput=boundingbox at PERMANENT binput=islands at PERMANENT
output=ocean operator=not

and export the resulting vector shape to a shapefile. Plotting ocean in
GRASS

      d.vect map=ocean at PERMANENT type=area fcolor=indigo

yields the image shown on
http://www.flickr.com/photos/39340925@N07/3616682800 (which looks exactly
like what I was looking for).

However after importing it into R with rgdal and plotting
> ocean <- readOGR("myfolder", "myfilename")
> plot(ocean, col="lightblue")

I get an image where the entire background is blue - not just the ocean,
but the islands as well.

Looking at the SpatialPolygonsDataFrame, I see that the polygons are not
marked as holes. But changing this does not seem to make any difference:

> boundary <- which.max(sapply(ocean at polygons, function(x)
x at Polygons[[1]]@area))
> for (i in ((1:length(ocean at polygons))[-boundary]))
ocean at polygons[[i]]@Polygons[[1]]@hole <- T

It appears as if the non-hole polygon is drawn fully filled, irrespective
of any holes, and then the 'holes' are drawn on top. How might i get around
that?

Thanks in advance for any suggestions,
Bjarke Christensen



More information about the R-sig-Geo mailing list