Change order of drawing geom_polygon, ggplot

Question

I am using R / ggplot to create a global map of the world with the world's earthquake zones. I am using a modified shapefile from here: http://gmo.gfz-potsdam.de/.

There are 9 layers in the plot and over 4400 groups. However, some of the layers should have holes, but appear not to. So to get around this I need to render these polygons first, such that they appear up the back. How do I change the render order of the polygons when using geom_polygon?

require("rgdal") # requires sp, will use proj.4 if installed
require("maptools")
require("ggplot2")
gpclibPermit() # required for fortify method
require("plyr")

fill_colours <- c("#dddddd", "#8FC463FF", "#46A742FF", "#FFFC00FF", "#F2B212FF", "#E9726AFF", "#D92C50FF", "#FF1B23FF", "#B1681FFF")

eq.shp = readOGR(dsn=".", "PGA_PolygonsPY v2_region")

# Create a discrete scale, there are probably more efficient ways of doing this
eq.shp$MAXVALUE[eq.shp$MAXVALUE==10.0] = 9
eq.shp$MAXVALUE[eq.shp$MAXVALUE== 4.8] = 8
eq.shp$MAXVALUE[eq.shp$MAXVALUE== 4.0] = 7
eq.shp$MAXVALUE[eq.shp$MAXVALUE== 3.2] = 6
eq.shp$MAXVALUE[eq.shp$MAXVALUE== 2.4] = 5
eq.shp$MAXVALUE[eq.shp$MAXVALUE== 1.6] = 4
eq.shp$MAXVALUE[eq.shp$MAXVALUE== 0.8] = 3
eq.shp$MAXVALUE[eq.shp$MAXVALUE== 0.4] = 2
eq.shp$MAXVALUE[eq.shp$MAXVALUE== 0.2] = 1
eq.shp$MAXVALUE  <- sprintf("%d", eq.shp$MAXVALUE)

eq.shp@data$id = rownames(eq.shp@data)
eq.points = fortify(eq.shp, region='id')

# Join the data frames back together
eq.df = join(eq.points, eq.shp@data, by='id')

p <- ggplot()
p <- p + geom_polygon(data = eq.df, aes(long, lat, group=group, fill=MAXVALUE))
p <- p + scale_fill_manual(values=fill_colours) 
p <- p + coord_equal(ratio = 1, xlim=c(65, 85), ylim=c(35, 45))
p    

This code produces the following plot.

If you compare the above plot with the original data data (plotted using GIS software), you can see the following artifacts: - The area around point (Long=77, Lat=42) in the original plot is not present. This is because the brown polygon is rendered last. - Likewise for the large green area around (Long=80, Lat=38), it is all yellow. This is because the yellow polygon is rendered after the green polygon. - I am not too concerned about the missing holes, e.g. area near (Long=74, Lat=39).

Answer 1

Okay, I have been able to solve this. Although it's a bit of a workaround and the render time is increased, it does meet my requirements.

Instead of trying to sort / reorder the rendering of the polygons using ggplot, I have forced this by using a for loop and to render the layers as I need them.

The solution:

Replace the following code:

p <- p + geom_polygon(data = eq.df, aes(long, lat, group=group, fill=MAXVALUE))

with the following:

for (i in seq(9, 1, -1)) {
    p <- p + geom_polygon(data = eq.df[eq.df$MAXVALUE==sprintf("%d", i),], aes(long, lat, group=group, fill=MAXVALUE))
}

This causes the polygons to be rendered in the order I wished it to. Then I get something very similar to the original:

However, I can see that I have introduced new issues, there are some polygons that don't have holes, so rendering them has causes them to cover some of the other layers. So I will need to create holes in the polygons to render them correctly.

Thank you to everyone who left comments.

Update: The solution above works to reorder the rendering of the polygons, but it does not entirely solve the underlying problem. It appears that you would need a more sophisticated algorithm to do so. I also tried ordering by the area of the polygon and then weighted area of the polygon, but neither worked to achieve and ideal solution. The final solution would require some manual intervention.