0

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

polygons plotted

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?

ggplot with extra polygons

1

2 Answers 2

2

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)

enter image description here

1
  • Yeah, I need to switch over to sf more generally. This project is convincing me of that.
    – jebyrnes
    Mar 2, 2018 at 23:05
0

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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.