[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