Shiny app that generates leaflet map out of a shapefile and CSV too slow; crashes with large shapefiles

I've been trying to build an app that lets me explore origin-destination flow data. It takes in a csv that has origin zone numbers, destination zone numbers, and a bunch of fields of information, plus a shapefile of the zones as polygons and lets me depict the flows to or from a specific zone. I've actually got a Shiny app that works, but it's either far too inefficient or maybe this is not quite the right tool for this job. I'll paste the code below.

My problem is that it seems to be too memory-intensive. For a small number of zones it works but it's extremely slow (takes ages to load initially and to update each time one of the inputs is changed). The bigger issue is that I normally work with large zone systems, with hundreds or even thousands of zones, with millions total vertices. When I try to apply this to one of those large zone systems it does load initially but after two or three movements (zooming or panning) or just one input change RStudio crashes.

So I was wondering if there's any way to make this more efficient or if RStudio/Shiny are simply not meant to be used for something like this. Keep in mind that I'm pretty much clueless about optimization so I'd appreciate it if you dumbed down the explanation a bit.

This is the app that works:

library(shiny)
library(rgdal)
library(sp)
library(leaflet)
library(RColorBrewer)

shape <- readOGR(dsn=".", layer="zones")
shape <- spTransform(shape, CRS("+init=EPSG:4326"))

Reading the shapefile, then the data. 'i' and 'j' are origin and destination zone numbers. I also need the list of possible fields to map, and it has to be made up of fields that are entirely numerical. I haven't figured out how to subset that without the loop. Any help on that would be appreciated.

values <- read.csv("values.csv")
zonelist <- sort(unique(c(values$i, values$j)))
zones <- length(zonelist)
fieldlist <- as.vector(character())
fields <- 0
for (valname in colnames(values)){
  if (valname != "i" &
      valname != "j" &
      is.numeric(values[[valname]])
      ) {
        fields <- fields + 1
        fieldlist[fields] <- valname
        }
}

And I'm putting them in a matrix mm. Again, any more efficient way to do this would be welcome.

mm <- array(numeric(), dim=c(fields,zones,zones))
mm[,,] <- 0
valarray <- array(numeric(), dim=fields*zones*zones)
valarray[] <- 0
for (k in 1:length(values)){
  if (colnames(values)[k] %in% fieldlist) {
    valarray[(values$j - 1) * zones * fields
             + (values$i - 1) * fields
             + which(colnames(values)[k] == fieldlist)
             ] <- values[[k]]
  } 
}
mm[,,] <- valarray

The interface is simple: choose from/to the reference zone, choose hte reference zone, choose the numerical field to map. Show a leaflet map of the flows.

ui <- fluidPage(
   titlePanel("Interzonal values mapping"),
   sidebarLayout(
      sidebarPanel(
         selectInput(inputId = "ft", label = "From/to zone", choices = c("from", "to")),
         selectInput(inputId = "zn", label = "Reference zone", choices = zonelist),
         selectInput(inputId = "fl", label = "Field to map", choices = fieldlist),
         width = 2
      ),
   mainPanel(leafletOutput("mapped_zonelist", height = 800),
             width = 10)
   )
)

And the server picks up the flows to map as a vector (vs) from the matrix mm depending on the inputs and formats the zone map according to the flow values.

server <- function(input, output) {
   m <- reactive({
     if (input$ft == "from"){
        vs <- mm[which(input$fl == fieldlist), as.numeric(input$zn), ]
       } else {
        vs <- mm[which(input$fl == fieldlist), , as.numeric(input$zn)]
       }
     
     c_weight <- numeric()
     c_weight[1:zones] <- 1
     c_weight[as.numeric(input$zn)] <- 2
     c_color <- character()
     c_color[1:zones] <- "blue"
     c_color[as.numeric(input$zn)] <- "black"

     leaflet() %>%
       addTiles() %>%
       addPolygons(data=shape,
                   fillOpacity = 0.7,
                   opacity = 0.8,
                   weight = c_weight[shape@data$n],
                   color = c_color[shape@data$n],
                   fillColor = colorNumeric(palette = brewer.pal(10, "RdBu"), domain = vs)(vs[shape@data$n]),
                   popup = paste0("Zone ", input$zn, " to zone ", as.character(shape@data$n),
                                  "<br>",
                                  input$fl, " = ", vs)
                  )
   })
   output$mapped_zonelist <- renderLeaflet(m())
}

shinyApp(ui = ui, server = server)

Solution

Depending on the size/complexity of your shapefile, the first thing I would do is simplify it. Load your shapefile as normal, then:

install.packages("rmapshaper")
shape_simp <- rmapshaper::ms_simplify(shape, keep = 0.05)

Experiment with the keep = argument. The greater the number, the greater the complexity of the resultant shapefile. This might not solve your problem entirely, but it's a start.