2

I've been struggling with this problem for a few days. I have a shapefile that contains zip code polygons. I'm trying to plot a certain portion of it with a transparent color (alpha=0.7), but some of the polygons at the edge of the plot end up with no color. If alpha=1.0, all polygons are fine. The shapefile is stored as a Large SpatialPolygonsDataFrame. 

# Code that you can use (Thanks jbaums)

library(rgdal);
download.file('http://www.naturalearthdata.com/http//www.naturalearthdata.com/do‌​wnload/110m/cultural/ne_110m_admin_0_countries.zip', {f <- tempfile()}); 
unzip(f, exdir=tempdir()); 
shp <- readOGR(tempdir(), 'ne_110m_admin_0_countries'); 
plot(shp, col='#FF000070', xlim=c(-100, 25))

For reference, my code looks like this :

#~~~ Read zip shapefile
  in_zip.shp  = "./Data/USZIP11_wSTCntyLLFIPS.shp"
  zip_poly = readShapePoly(in_zip.shp)

#~~~ Set-up the mapping parameters
  xlim = c(-82.0, -80.7)
  ylim = c(25.5, 26.8)

#~~~ Works well
  output_figure <- "./Good.png"
  png(output_figure, width=5, height=4, units="in", res=1200)
  par(mar=c(1.5, 1.8, 2.0, 0.7))
  plot(xlim, ylim, type="n", xlim=xlim, ylim=ylim, axes=FALSE, xlab="", ylab="")

  plot(zip_poly, col=rgb(1, 0, 0, alpha=1), xlim=xlim, ylim=ylim, lty=1.3, add=T)
  box()
  title(main="Good")
  dev.off()

#~~~ Doesn't work well
  output_figure <- "./Bad.png"
  png(output_figure, width=5, height=4, units="in", res=1200)
  par(mar=c(1.5, 1.8, 2.0, 0.7))
  plot(xlim, ylim, type="n", xlim=xlim, ylim=ylim, axes=FALSE, xlab="", ylab="")

  plot(zip_poly,col=rgb(1,0,0,alpha=0.7),xlim=xlim,ylim=ylim,lty=1.3,add=T)
  box()
  title(main="Bad")
  dev.off()

Any tips would be appreciated. I attached figures so you can see what problem I get with my code. Thanks!

good.png bad.png

jgadoury
  • 293
  • 2
  • 13
  • This behaviour can be replicated with: `library(rgdal); download.file('http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/cultural/ne_110m_admin_0_countries.zip', {f <- tempfile()}); unzip(f, exdir=tempdir()); shp <- readOGR(tempdir(), 'ne_110m_admin_0_countries'); plot(shp, col='#FF000070', xlim=c(-100, 25))` – jbaums Jun 03 '14 at 23:40
  • The issue is not specific to your data, so for generality and reproducibility, perhaps consider editing your question such that it uses the code I included in the comment above. It's up to you, but it might attract more feedback if it's a simple (writing to png, modifying `par`, etc is unnecessary clutter), reproducible example of your problem. I could then add plots to the question if you wanted. :) – jbaums Jun 03 '14 at 23:52
  • 1
    Can't reproduce it on Mac OS X 10.7 with R 2.14... – plannapus Jun 04 '14 at 07:36
  • 1
    Strange issue. I reproduced the problem with @jbaums code on Windows with R 2.15. The problems disappears if the xlim is changed to eg. c(-100,95). As a workaround, I clipped the map to the desired xlim and ylim and then worked alright. For clipping, see [this](http://stackoverflow.com/a/13986029/2516066) answer. – koekenbakker Jun 04 '14 at 14:29
  • For the record, I'm on Windows 8.1 x64 with R 3.1.0. Just tried it on Win 7 x32 with R 3.0.2 with the same results. RStudio on both. Also noticed that the url in my comment above mistakenly has white space in it (in case anyone else is trying to run the code). – jbaums Jun 04 '14 at 14:38
  • 1
    Thanks jbaums, I edited my question and added your code. I left mine there just for reference. – jgadoury Jun 04 '14 at 14:45
  • And for the record, I am also using R 3.0.2 with RStudio 0.98.501 in Windows Server 2008 R2 Standard 64-bit (VMware Virtual Platform) – jgadoury Jun 04 '14 at 14:50
  • @koekenbakker Thanks for this! I used clipping and it worked perfectly. I'll definitely put that in my code, but I'm still wondering why this issue happened in the first place.. – jgadoury Jun 04 '14 at 15:12

1 Answers1

1

First, I don't know the reason of the problem. The fact that it works for regular colours but not for colours with alpha channel, seems strange to me.

Second, here is an example of a work-around using gIntersection of the rgeos package to clip the polygons to the size of interest, i.e. the xlim and ylim. Plotting the clipped map seems to work for alpha colours. 

library(rgdal)
library(rgeos)  # for gIntersection
library(raster) # for extent()

# download world map
download.file('http://www.naturalearthdata.com/http//www.naturalearthdata.com/download/110m/cultural/ne_110m_admin_0_countries.zip', {f <- tempfile()}); 
unzip(f, exdir=tempdir()); 
shp <- readOGR(tempdir(), 'ne_110m_admin_0_countries'); 

# define map limits
xlim = c(0, 90)
ylim = c(40, 60)

# red with alpha
mycol = rgb(1,0,0,0.5)

# plot world map
plot(shp)

# create rectangle from xlim and ylim
CP <- as(extent(c(xlim,ylim)), "SpatialPolygons")
proj4string(CP) <- proj4string(shp)

# add to plot
plot(CP, add=T, col=rgb(1,0,0,0.5))

# set up window for 2 plots, comparing 2 methods
par(mfrow=c(2,1))

# plot map with xlim and ylim: does not work the way we want
plot(shp, xlim=xlim, ylim=ylim, col=mycol, axes=T, las=1)

# clip map to xlim and ylim, then plot: works
out <- gIntersection(shp, CP, byid=TRUE)
plot(out, col=mycol, axes=T, las=1)
koekenbakker
  • 3,524
  • 2
  • 21
  • 30
  • Thanks a lot for that solution, worked like a charm. I previously tried something that consisted of changing the coords of the polygons to never be out of `xlim` and `ylim`, but then I ran in the issue that some polygons became actually lines at the edge of the box and R didn't like that. – jgadoury Jun 04 '14 at 15:16