ggplot2 adding extra polygons

Question

So, I'm working with a dataset of polygons generated at different levels of sensitivity to error from a citizen science project. I have been trying to visualize how my sensitivity threshold changes results, but in plotting have noticed odd behaviour in ggplot2, inserting polygons where there should be none. The comparison between plotting the native SpatialPolygonsDataSet and the tidy generated data frame is particularly telling. What's going on here, and how can I fix it?

I start with reading in the data and doing some light processing

library(tidyverse)
library(broom)
library(sp)
library(spdplyr)
library(ggplot2)
library(rgdal)
library(raster)
library(rgeos)
library(lubridate)


#read in spatialPolygonsDataFrame (shapefile)
#from https://www.nceas.ucsb.edu/~byrnes/floating_forests/ff_polys_proj.rds
classifications <- readRDS("./ff_polys_proj.rds")

classifications <- classifications %>% 
  mutate(scene_timestamp = parse_date_time(scene_timestamp, orders="ymdHMS"),
         quarter = quarter(scene_timestamp, with_year=TRUE, fiscal_start=11))

(note, this all takes some time as the file is big-ish)

Here's a simple plot of the result focusing on one area filtering down to those polygons with a threshold of 6.

#plot the north coast at a threshold of 6
plot(classifications %>% filter(threshold==6), 
     xlim=c(-124.001673, -122.772273), 
     ylim=c(38.115699,39.176660) )

OK, now, here's a data frame instead

#convert to data frame
classifications_df <- classifications %>% tidy
classifications@data$id <- rownames(classifications@data)
classifications_df <- left_join(classifications_df, classifications@data)

which can then be ggplotted

##ggplot the north coast
ggplot() + 
  geom_polygon(data=classifications_df %>% filter(threshold==6),  
               mapping=aes(x=long, y=lat, group=group)) +
  theme_bw(base_size=14) +
  xlab("") + ylab("") + 
  xlim(-124.001673, -122.772273) + 
  ylim(38.115699,39.176660)

And now... why are there extra polygons where formerly there were none?

Answer 1

You have asked for a solution with SpatialPolygonsDataFrame, but plotting polygon data is now quite a bit easier using package sf and ggplot2::geom_sf. Try this:

# libraries required is considerably reduced. Note you need development version of ggplot2 (devtools::install_github('tidyverse/ggplot2'))
library(ggplot2)
library(lubridate)
library(sf)

#read in spatialPolygonsDataFrame (shapefile)
#from https://www.nceas.ucsb.edu/~byrnes/floating_forests/ff_polys_proj.rds
classifications <- readRDS('input/ff_polys_proj.rds')

# reduce columns and convert to sf object
classifications <- classifications[,c('scene_timestamp', 'threshold')] %>%
  st_as_sf()

# dplyr functions work with sf
classifications <- classifications %>% mutate(scene_timestamp = parse_date_time(scene_timestamp, orders = "ymdHMS"),
         quarter = quarter(scene_timestamp, with_year = TRUE, fiscal_start = 11))

# set limits within coord_sf, set axis labels to NULL rather than "", so that space is not reserved in the margins.
gp <- ggplot() +
  geom_sf(data = classifications %>% filter(threshold == 6), lwd = 0.3, fill = 'black') +
  coord_sf(xlim = c(-124.001673, -122.772273), ylim = c(38.115699,39.176660)) +
  theme(axis.title = NULL)

# save to file
ggsave('output/map-class.png', plot = gp)

Answer 2

OK, the answer appears to be to not rely on rownames as I have in the past, but rather get the IDs right from the polygons. This seems.... opaque. Is there a better solution that I'm missing here?

 classifications@data$id <- sapply(classifications@polygons, function(x) x@ID)