Change SP Object Polygon Rendering Order in R

The Problem

I have a geoJSON file that was made by combining many (as in ~250) geoJSON files each containing a single polygon in R.  The polygons are of varying size and can overlap in space.  When I open the file in a GIS, often the smaller polygons are hidden under the larger ones, making displaying the data challenging, if not impossible.

Map showing few polygons when there should be many

I know there are more polygons than are visible in this area, so they must be hiding under the larger polygons.

Map showing larger polygons obscuring smaller polygons

With transparency set to 60% for each polygon (see the Draw Effects dialog in QGIS for the fill symbol layer), you can see that smaller polygons are hiding under larger ones.

The Goal

I would prefer that the polygons stack up so that the largest is on the bottom and the smallest is on the top.  This requires that I change the rendering order based on the area of each polygon.

The Quick Fix

QGIS has an option to control the rendering order.  Open the Layer Properties; go to the Style tab; check the box for “Control feature rendering order”; click the button on the right with the A Z and an arrow to enter the expression you want (I would order by area for example).

Why isn’t this the right solution for me?  I’m sharing a user-contributed dataset.  One of the goals is that anyone can use it.  When polygons are obscured, it makes the dataset just a little harder to use and understand, which means people won’t like using it.  Another goal is that anyone with a reasonable understanding of GIS can contribute.  If I have to write a bunch of instructions on how to visualize the data before they can add to the dataset, people are less likely to contribute.

Map showing all the polygons expected.

Now I can see all of the polygons because the larger ones are on the bottom and the smaller ones are on top.

My Solution

Hat tip to Alex Mandel and Ani Ghosh for spending a couple of hours with me hammering out a solution.  We used R because I already had a script that takes all of the polygons made by contributors and combines them into one file.  It made sense in this case to add a few lines to this post-processing code to re-order the polygons before sending the results to the GitHub repository.

What you need to know about rendering order & SP Objects

The order in which items in an SP object are rendered is controlled by the object’s ID value.  The ID value is hidden in the ID slot nested inside the polygons slot.  If you change these values, you change the order items are rendered.  ID = 1 goes first, ID =2 goes on top of 1, 3 goes on top of 2, and so on.  So for my case, assigning the IDs based on the area will get me the solution.

How

This Stack Exchange Post on re-ording spatial data was a big help in the process.  Note that every SP object should have the slots and general structure I used here.  There’s nothing special about this dataset.  If you’d like the dataset and another copy of the R code, however, it is in the UC Davis Library’s AVA repository.

#load the libraries you'll need
library(raster)
library(geojsonio)
library(rgdal)

### FIRST: Some context about how I made my dataset in the first place

# search in my working directory for the files inside the folders 
# called "avas" and "tbd"
avas <- list.files(path="./avas", pattern = "*json$", full.names = "TRUE")
tbd <- list.files(path="./tbd", pattern = "*json$", full.names = "TRUE")

#put the two lists into one list
gj <- c(avas, tbd)

#read all the geojson files & create an SP object
vects <- lapply(gj, geojson_read, what="sp")

#combine all the vectors together. bind() is from the raster package.
#probably could just rbind geojson lists too, but thats harder to plot
all <- do.call(bind, vects)

#Change any "N/A" data to nulls
all@data[all@data=="N/A"]<- NA


### SECOND: How I did the sorting

#Calculate area of polygons - needed for sorting purposes
# the function returns the value in the area slot of each row
all@data$area<-sapply(slot(all, "polygons"), function(i){slot(i, "area")})

#add the row names in a column - needed for sorting purposes
all$rows<-row.names(all)

#Order by area - row names & area are needed here
# decreasing = TRUE means we list the large polygons first
all<-all[match(all[order(all$area, decreasing = TRUE),]$rows, row.names(all@data)),]

#add new ID column - essentially you are numbering the rows 
# from 1 to the number of rows you have but in the order of 
# largest to smallest area
all$newid<-1:nrow(all)

#assign the new id to the ID field of each polygon
for (i in 1:nrow(all@data)){
 all@polygons[[i]]@ID<-as.character(all$newid[i])}

#drop the 3 columns we added for sorting purposes (optional)
all@data<-all@data[,1:(ncol(all@data)-3)]

#write the data to a file in my working directory
geojson_write(all, file="avas.geojson", overwrite=TRUE, convert_wgs84 = TRUE)
Advertisements

About micheletobias

I lead two lives - one as an artist and the other as a scientist. More and more I'm finding my two worlds colliding, and it's not the disaster you might expect. View all posts by micheletobias

I'd love to hear what you think. Leave a comment!

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: