3

I am plotting multiple shapefiles using spplot. Here's a data to construct that

library(raster)
library(randomcoloR)

my.shp <- getData('GADM', country = 'BRA', level = 2)
my.shp$ID<- 1:nrow(my.shp)

My data consists of a variable X for 10 years as shown where each column is a year

df <- matrix(sample(100:5000, 55040, replace = T), nrow = 5504, ncol = 10)
df <- data.frame(ID = 1:nrow(my.shp), df)

my.dat <- merge(my.shp, df, by = "ID")

variable.names <- paste0("X",1:10)

spplot(my.dat, rev(variable.names), col = NA, at = seq(from = 100, to = 5000, by = 500), 
          col.regions = distinctColorPalette(length(seq(from = 100, to = 5000, by = 500))),
          main = list(label = "TEST")) 

My problem is this plot takes so much time (around an hour) to get plotted and was wondering if there is something inherently wrong in the code itself that it is taking too long to plot. My laptop has a 32 GB RAM.

Thanks

89_Simple
  • 3,393
  • 3
  • 39
  • 94

2 Answers2

2

I haven't compared this plot to your spplot because I don't want to spend an hour waiting for it.

Instead I'm proposing to use library(mapdeck) to plot an interactive map, which takes a matter of seconds.

Two things to note

  1. You need a Mapbox Access token
  2. You need to convert the sp object to sf
library(raster)

my.shp <- getData('GADM', country = 'BRA', level = 2)
my.shp$ID <- 1:nrow(my.shp)

df <- matrix(sample(100:5000, 55040, replace = T), nrow = 5504, ncol = 10)
df <- data.frame(ID = 1:nrow(my.shp), df)

my.dat <- merge(my.shp, df, by = "ID")


library(sf)
sf <- sf::st_as_sf( my.dat )

library(mapdeck)

set_token( "YOUR_MAPBOX_TOKEN" )

mapdeck() %>% 
  add_sf(
    data = sf
    , fill_colour = "GID_2"
    )

enter image description here

SymbolixAU
  • 25,502
  • 4
  • 67
  • 139
  • Thanks though it looks like a good alternative, I cannot unfortunately award my bounty since I am interested in the solution to my exact problem. I will do use your approach for other stuff. Thanks for your time and effort. – 89_Simple Jan 23 '19 at 10:50
  • Interestingly, I see a blank page using this approach. Any reason why that would be the cause? I tried using both the public as well as private token – 89_Simple Jan 23 '19 at 10:52
  • Try opening it in a browser. See here too https://symbolixau.github.io/mapdeck/articles/issues.html – SymbolixAU Jan 23 '19 at 10:58
  • Okay thanks. Looks like I can only open rstudio on a web browser in linux. I only know how to work on windows. Thanks anyway – 89_Simple Jan 23 '19 at 11:06
  • The browser should open no matter which platform you're on. – SymbolixAU Jan 23 '19 at 20:33
  • "RStudio is available in open source and commercial editions and runs on the desktop (Windows, Mac, and Linux) or in a browser connected to RStudio Server or RStudio Server Pro (Debian/Ubuntu, RedHat/CentOS, and SUSE Linux)."..so does it mean R studio server is not available for windows? – 89_Simple Jan 24 '19 at 13:32
  • 1
    I don't mean open RStudio in a browser, I mean open the plot in a browser. See https://stackoverflow.com/a/52470812/5977215 – SymbolixAU Jan 24 '19 at 19:25
1

Are you willing/able to switch to sf instead of sp?

The sf plot function is considerably faster than spplot, although the layout differs a bit.

library(sf)
my.dat_sf <- st_as_sf(my.dat)
plot(my.dat_sf[rev(variable.names)], max.plot=10, breaks=c(seq(from = 100, to = 5000, by = 500),5000),
     pal = distinctColorPalette(length(seq(from = 100, to = 5000, by = 500))),
     main = "TEST", border=NA, key.pos=4)

Additionally, you could try to simplify the polygon with rmapshaper::ms_simplify() for Spatial*-objects or sf::st_simplify() for SimpleFeatures, which lets you reduce the object size by quite a bit, depending on the given dTolerance. Thus plotting, will also be faster with simplified polygons.

The original SpatialPolygon:

format(object.size(my.dat_sf), units="Kb")

"25599.2 Kb"

and a simplified SimpleFeature:

dat_sf_simple <- st_transform(my.dat_sf, crs = 3035)
dat_sf_simple <- st_simplify(dat_sf_simple, dTolerance = 1000, preserveTopology = T)
dat_sf_simple <- st_transform(dat_sf_simple, crs = 4326)
format(object.size(dat_sf_simple), units="Kb")

"7864.2 Kb"

Plot the simplified SimpleFeature, which takes about 1 minute on my machine with 8GB RAM.

plot(dat_sf_simple[rev(variable.names)], max.plot=10, breaks=c(seq(from = 100, to = 5000, by = 500),5000),
     pal = distinctColorPalette(length(seq(from = 100, to = 5000, by = 500))),
     main = "TEST", border=NA, key.pos=4)

You could also try out with ggplot2, but I am pretty sure the most performant solution will be the sf plot.

library(ggplot2)
library(dplyr)
library(tidyr)

dat_sf_simple_gg <- dat_sf_simple %>% 
  dplyr::select(rev(variable.names), geometry) %>% 
  gather(VAR, SID, -geometry)

ggplot() +
  geom_sf(data = dat_sf_simple_gg, aes(fill=SID)) + 
  facet_wrap(~VAR, ncol = 2) 
SeGa
  • 9,454
  • 3
  • 31
  • 70
  • Okay. I have two follow up questions: my solution allows me to have a common legend across all the maps. How do I insert a common legend in your solution? Second If I want to add the state boundary on top, in my solution I can just do `latticeExtra::layer(sp.polygons(my.boundary, lwd = 0.05))`. How does that work in your case? – 89_Simple Jan 23 '19 at 12:49
  • To see the state boundaries, remove `border=NA`. And to add a legend, add `key.pos=4` in the plot call. – SeGa Jan 23 '19 at 13:11
  • Thanks. What I meant was if I want to add another polygon on top of this map, how do I do it. Doing border = NA will only remove the current boundary of the municipalities on this – 89_Simple Jan 23 '19 at 14:00
  • Could you include the call to `latticeExtra::layer` in your example. I am not able to make that plot and therefore I dont really know hwo the result should look like. But sf plot also takes the argument `add`. But I dont really know how it will add to the plot, if the plot has "facets" when plotting multiple columns. – SeGa Jan 23 '19 at 16:04