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

[Roger Bivand]

On Thu, 11 Jun 2009, Bjarke Christensen wrote:

>
> 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.

As you observe, the approach you are taking does not work, as the R 
graphics devices work by over-painting in layers. To see the image, it has 
to be painted after the enclosing rectangle. Holes are painted by 
re-painting the hole in a chosen background colour, which by default is 
"transparent", so you just get lots of blue.

Why not do something like

o <- overlay(myDEM, islands)

(untried) to get just the raster cells within the island polygons as a 
SpatialPixelsDataFrame object and image() that, possibly setting the 
background to a suitable colour (NAs will get transparent by default). 
This works with the graphics system, rather than trying to work round it - 
it doesn't "remember" that there is anything on the canvas that should be 
protected from overpainting, so it is safer just to paint what needs 
painting. I can also imagine painting first with reduced opacity (or 
intensity) in a different palette, then overpainting just the islands with 
full intensity in the target palette, which might be more visually 
pleasing than just flat blue sea.

Hope this helps,

Roger

>
> 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
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no