R/leaflet_plot.R

In NERC-CEH/ukem: Package Title (one line, title case)

## some code to make a plot to check outputs in "./data/output/"

packs <- c("terra","stringr","plyr","ggplot2","data.table","stats","readr","sf","readxl","targets","tidyverse","ggforce","spCEH","leaflet","shiny","shinydashboard","leaflet.extras", "viridis")
lapply(packs, require, character.only = TRUE)

source("./R/ukem.R")

## set up sector names

# shiny
ui <- bootstrapPage(
  tags$style(type = "text/css", "html, body {width:100%;height:100%}"),
  leafletOutput("rasmap", width = "100%", height = "100%"),
  absolutePanel(class = "panel panel-default", fixed = TRUE,draggable = TRUE,
                top = 20, left = 100,
                style="padding-left: 5px; padding-right: 5px; padding-top: 5px; padding-bottom: 5px",
                width = "auto", height = "auto",
                style = "opacity: 0.75",
                selectInput("year", "Year: ", c(2010:2019), selected = 2019),
                selectInput("ghg", "GHG: ", c("ch4","co2","n2o"), selected = "ch4"),
                selectInput("sector", "SNAP: ", c(1:11,"Total"), selected = "Total"),
                sliderInput("opacity", "Opacity %:",min=1,max=100,value=90),
                plotOutput("emisKt", width="100%",height = "100%"),
                plotOutput("emisRel", width="100%",height = "100%")

  )
)


server <- function(input, output,session) {

  # input totals by SNAP
  l_ukem_input <- eventReactive(c(input$year, input$ghg), {

    dt_ts_uk <- tar_read(dt_ts_uk)
    dt_ts_ie <- tar_read(dt_ts_ie)
    dt_ts <- rbindlist(list(dt_ts_uk, dt_ts_ie), use.names = T, fill=T)
    dt_ts <- dt_ts[pollutant == input$ghg & year == input$year]
    dt <- dt_ts[, .(input_emis = sum(emission,na.rm=T)), by = .(snap_id, pollutant, year)]
    dt[,input_emis := units::set_units(input_emis, Gg / yr)]

    list(dt = dt)
  })

  # output totals per SNAP
  l_ukem_output <- eventReactive(c(input$year, input$ghg), {

    molWt <- unit_conversion(pollutant=input$ghg, "mol")$molWt

    fname <- paste0("./data/output/s_F_",input$ghg,"_",input$year,".tif")

    s <- rast(fname)
    # cell areas
    r_area <- terra::cellSize(s[[1]], unit="m")

    dt <- data.table(snap_id=1:11, output_emis = unlist(global(s[[1:11]]*r_area, sum, na.rm=T) * molWt * 60 * 60 * 8760 / 1000000000))
    dt[,output_emis := units::set_units(output_emis, Gg / yr)]

    list(dt = dt, s = s, r_area = r_area, molWt = molWt)
  })

  # output raster, converted to tonnes per cell yr-1
  r_ukem_output <- reactive({

    # subset
    if(input$sector != "Total"){

      r <- l_ukem_output()$s[[as.numeric(input$sector)]]

    }else{

      r <- app(l_ukem_output()$s[[1:11]], fun="sum", na.rm=T)

    }

    # convert from mol m-2 s-1 to t yr-1
    r_t_s <- (r * l_ukem_output()$r_area * l_ukem_output()$molWt * 60 * 60 * 8760)/1000000
    LL_total <- unlist(global(r_t_s, sum, na.rm=T))

    # reproject for plot, rescale back to LL total
    r_t_s <- suppressWarnings(terra::project(r_t_s, "EPSG:3857"))
    r_t_s <- (r_t_s / unlist(global(r_t_s, sum, na.rm=T))) * LL_total

    r_t_s <- raster(r_t_s)

    if(cellStats(r_t_s,sum) != 0) r_t_s[r_t_s==0] <- NA

    list(r_t_s = r_t_s)

  })

  # base map
  output$rasmap <- renderLeaflet({
    leaflet() %>%
      addProviderTiles("CartoDB.Voyager", layerId = "rastile",
                       options = providerTileOptions(minZoom = 4, maxZoom = 12))%>%
      setView(lng = -2, lat= 54, zoom = 6)
  })

  # add map depending on SNAP
  observe({

    pal <- colorNumeric(viridis(8), domain = values(r_ukem_output()$r_t_s) ,na.color = "transparent")

    proxy <- leafletProxy("rasmap")
    proxy %>% removeTiles(layerId="rasimg") %>%
      removeControl(layerId="rasleg") %>%
      addRasterImage(r_ukem_output()$r_t_s, colors = pal, project=F, layerId="rasimg", opacity = (input$opacity)/100) %>%
      clearPopups()%>%
      addLegend(pal = pal,
                values = values(r_ukem_output()$r_t_s),
                position = "topleft",
                title = "t yr-1", layerId = "rasleg")
  })


  # plot for total emissions
  output$emisKt <- renderPlot({

    ggplot(data=l_ukem_output()$dt, aes(x=snap_id, y=output_emis))+
      geom_bar(stat="identity")+
      ylab("Emission (Kt)")+
      theme_bw()

  }, width=300, height=250)

  # plot for inputs vs outputs (to 1)
  output$emisRel <- renderPlot({

    dt1 <- copy(l_ukem_input()$dt)
    dt2 <- copy(l_ukem_output()$dt)
    x <- dt1[dt2, on = "snap_id"]
    x[, acc_for := units::drop_units(output_emis / input_emis)]

    ggplot(data=x, aes(x=snap_id, y=acc_for))+
      geom_bar(stat="identity")+
      ylab("Emission: Output / Input")+
      geom_hline(yintercept=1, linetype="dashed", colour="red")+
      theme_bw()

  }, width=300, height=250)


  # pop up click event
  observeEvent(input$rasmap_click, {

    click <- input$rasmap_click

    xy <- st_as_sf(data.table(lng = click$lng, lat = click$lat), coords = c("lng","lat"), remove = FALSE)
    st_crs(xy) <- "EPSG:4326"
    xy <- st_transform(xy, "EPSG:3857")
    rasval <- extract(r_ukem_output()$r_t_s, xy)
    content <- paste("<strong>","emission =", round(rasval, 3)," t","</strong>", "<br>",
                     "Lon =", round(click$lng, 2),"<br>",
                     "Lat =", round(click$lat, 2))

    if (is.na(rasval))
      return()

    proxy <- leafletProxy("rasmap")
    proxy %>% clearPopups() %>%
      addPopups(click$lng, click$lat, content)
  })

}

runApp(shinyApp(ui, server), launch.browser = TRUE)

#####################################################################
#####################################################################