1. Adds map plot of model output to priors app

CHANGELOG.md

@@ -15,6 +15,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Clearer documentation on types expected by the `mesh_builder` tool - [PR #101](https://github.com/4DModeller/fdmr/pull/101)
 - Checks on the data types being passed into the `mesh_builder` tool - [PR #101](https://github.com/4DModeller/fdmr/pull/101)
 - Ability to plot either polygon or point data on Leaflet map of `mesh_builder` tool - [PR #101](https://github.com/4DModeller/fdmr/pull/101)
+- The ability to plot model predictions on a `leaflet` map in the our [Interactive priors Shiny app](https://4dmodeller.github.io/fdmr/articles/priors_app.html) - [PR #147](https://github.com/4DModeller/fdmr/pull/147)
 
 ### Fixed
 

DESCRIPTION

@@ -54,7 +54,8 @@ Imports:
     shinybusy,
     curl,
     promises,
-    fmesher
+    fmesher,
+    purrr
 Suggests: 
     bookdown,
     knitr,
@@ -66,4 +67,3 @@ Suggests:
     rcmdcheck
 Config/testthat/edition: 3
 VignetteBuilder: knitr
-

NAMESPACE

@@ -9,6 +9,7 @@ export(load_tutorial_data)
 export(mesh_builder)
 export(mesh_checker)
 export(mesh_to_spatial)
+export(model_viewer)
 export(numbers_only)
 export(parse_model_output)
 export(plot_barchart)

R/model_parse.R

@@ -43,3 +43,26 @@ parse_model_output <- function(model_output, measurement_data, model_type = "inl
         return(parse_model_output_bru(model_output = model_output, measurement_data = measurement_data))
     }
 }
+
+
+#' Create a prediction field from the parsed model output and the mesh
+#'
+#' @param var_a 
+#' @param var_b 
+#' @param mesh 
+#'
+#' @return data.frame
+#' @export
+create_prediction_field <- function(var_a, var_b, mesh) {
+    mod_proj <- fmesher::fm_evaluator(mesh)
+    xy_grid <- base::expand.grid(mod_proj$x, mod_proj$y)
+    A_proj <- INLA::inla.spde.make.A(mesh = mesh, loc = as.matrix(xy_grid))
+
+    z <- base::exp(base::as.numeric(A_proj %*%  var_a[1:mesh$n]) + base::sum(var_b))
+    base::data.frame(x = xy_grid[, 1], y = xy_grid[, 2], z = z)
+}
+
+
+create_raster <- function(dataframe, crs) {
+    raster::rasterFromXYZ(dataframe, crs = crs)
+}

R/shiny_priors.R

@@ -157,9 +157,16 @@ priors_shiny <- function(spatial_data,
                         shiny::selectInput(inputId = "plot_type", label = "Plot type:", choices = plot_choices, selected = plot_choices[1]),
                         shiny::plotOutput(outputId = "plot_model_out")
                     ),
+                    shiny::tabPanel(
+                        "Map",
+                        shiny::selectInput(inputId = "select_run_map", label = "Select run:", choices = c()),
+                        shiny::selectInput(inputId = "map_var_a", label = "Variable a:", choices = c()),
+                        shiny::selectInput(inputId = "map_var_b", label = "Variable b:", choices = c()),
+                        leaflet::leafletOutput(outputId = "map_out")
+                    ),
                     shiny::tabPanel(
                         "Code",
-                        shiny::selectInput(inputId = "select_run", label = "Select run:", choices = c()),
+                        shiny::selectInput(inputId = "select_run_code", label = "Select run:", choices = c()),
                         shiny::verbatimTextOutput(outputId = "code_out")
                     ),
                     shiny::tabPanel(
@@ -200,8 +207,32 @@ priors_shiny <- function(spatial_data,
             names(model_vals$run_params)
         })
 
+        model_summary_variables <- shiny::reactive({
+            names(model_vals$model_outputs)
+        })
+
         shiny::observe({
-            shiny::updateSelectInput(session = session, inputId = "select_run", choices = run_names())
+            shiny::updateSelectInput(session = session, inputId = "select_run_map", choices = run_names())
+            shiny::updateSelectInput(session = session, inputId = "select_run_code", choices = run_names())
+
+            if (!is.null(run_names())) {
+                # They'll all have the same variable names so we can just take the first
+                run_name <- run_names()[1]
+                var_names <- names(model_vals$parsed_outputs[[run_name]])
+
+                selected_a <- NULL
+                selected_b <- NULL
+                if ("mean_post" %in% var_names) {
+                    selected_a <- "mean_post"
+                }
+
+                if ("fixed_mean" %in% var_names) {
+                    selected_b <- "fixed_mean"
+                }
+
+                shiny::updateSelectInput(session = session, inputId = "map_var_a", choices = var_names, selected = selected_a)
+                shiny::updateSelectInput(session = session, inputId = "map_var_b", choices = var_names, selected = selected_b)
+            }
         })
 
         shiny::observeEvent(input$features, {
@@ -300,8 +331,10 @@ priors_shiny <- function(spatial_data,
                     function(model_output) {
                         # Run the model
                         run_no(run_no() + 1)
-                        model_vals$model_outputs[[run_no()]] <- model_output
-                        model_vals$parsed_outputs[[run_no()]] <- parse_model_output(
+                        run_label <- paste0("Run-", run_no())
+
+                        model_vals$model_outputs[[run_label]] <- model_output
+                        model_vals$parsed_outputs[[run_label]] <- parse_model_output(
                             model_output = model_output,
                             measurement_data = measurement_data
                         )
@@ -316,7 +349,6 @@ priors_shiny <- function(spatial_data,
                             "pg_ar1" = input$pg_ar1
                         )
 
-                        run_label <- paste0("Run-", run_no())
                         model_vals$run_params[[run_label]] <- run_params
 
                         if (write_logs) {
@@ -355,6 +387,27 @@ priors_shiny <- function(spatial_data,
             rownames = TRUE
         )
 
+        map_raster <- shiny::reactive({
+            if (length(model_vals$parsed_outputs) == 0) {
+                return()
+            }
+
+            data <- model_vals$parsed_outputs[[input$select_run_map]]
+            pred_field <- create_prediction_field(var_a = data[[input$map_var_a]], var_b = data[[input$map_var_b]], mesh = mesh)
+            create_raster(dataframe = pred_field, crs = sp::proj4string(spatial_data))
+        })
+
+        output$map_out <- leaflet::renderLeaflet({
+            if (is.null(map_raster())) {
+                return()
+            }
+
+            m <- leaflet::leaflet()
+            m <- leaflet::addTiles(m, group = "OSM")
+            m <- leaflet::addRasterImage(m, map_raster(), opacity = 0.9, group = "Raster")
+            m
+        })
+
         model_plot <- shiny::eventReactive(input$plot_type, ignoreNULL = FALSE, {
             if (length(model_vals$parsed_outputs) == 0) {
                 return()

vignettes/priors_app.Rmd

@@ -41,13 +41,13 @@ utils::head(covid19_data)
 
 # Create the mesh
 
-We can have a look at the mesh and change it interactively using the `mesh_builder` Shiny app. Using the app I came up with a mesh with the following parameters, we'll use this
-and pass it into the Setting Priors Shiny app below.
+We can have a look at the mesh and change it interactively using the `fdmr::mesh_builder` Shiny app. Using the app I came up with a mesh with the following parameters, we'll use this and pass it into the "Setting Priors Shiny app" below.
 
 ```{r createmesh, error=TRUE}
 initial_range <- diff(base::range(sp_data@data[, "LONG"])) / 3
 max_edge <- initial_range / 2
-mesh <- INLA::inla.mesh.2d(
+
+mesh <- fmesher::fm_mesh_2d_inla(
   loc = sp_data@data[, c("LONG", "LAT")],
   max.edge = c(1, 2) * max_edge,
   offset = c(initial_range, initial_range),
@@ -64,7 +64,7 @@ n_groups <- length(unique(covid19_data$week))
 
 # Set coordinates on data
 
-We will use the function `bru()` of package `inlabru` to fit the model. `bru` expects the coordinates of the data, thus we transform `covid19_data_bris` data set to a SpatialPointsDataFrame using the function `coordinates()` of the `sp` package.
+We will use the function [`inlabru::bru`](https://inlabru-org.github.io/inlabru/reference/bru.html) function to fit the model which expects the coordinates of the data, thus we transform `covid19_data` data set to a [`SpatialPointsDataFrame`](https://www.rdocumentation.org/packages/sp/versions/2.0-0/topics/SpatialPointsDataFrame-class) using the [`sp::coordinates`](https://www.rdocumentation.org/packages/sp/versions/2.0-0/topics/coordinates) function.
 
 ```{r, error=TRUE, eval=FALSE}
 sp::coordinates(covid19_data) <- c("LONG", "LAT")
@@ -74,10 +74,29 @@ sp::coordinates(covid19_data) <- c("LONG", "LAT")
 
 Now we have the filtered data we are ready to pass in the spatial 
 
-First start by selecting the variable to model. In this example we'll select `cases`, then select the features to add to the formula. At the bottom of the window you'll see the formula being constructed.
-Click on the checkboxes to select your formula and then click the `Run Model` button to run the model see it's output plot on the right of the window.
+## Selecting priors
+
+First start by selecting the variable to model. In this example we'll select `cases`, then select the features to add to the formula. At the bottom of the window you'll see the formula being constructed. Once you've setup the formula to your specifications click on the `Model` tab and click `Run`. Once the model run has finished you'll see a table of summarised model outputs.
+
+## Plotting model outputs
+
+Once the model has run you view plots of model outputs by clicking on the `Plot` tab. Each time you run a model with different sets or priors the parsed model output saved and can be compared with previous runs. 
+
+## Plotting model ouputs on a map
+
+To plot model predictions on a map click on the `Map` tab. This will use the `fdmr::create_prediction_field` and `fdmr::create_raster` functions to first create a `data.frame` of model predictions and then create a raster image from this data.
+
+## Saving model outputs
+
 By default the app will write out parameter sets and logs to the `fdmr/logs` directory in your home folder. If you want these logs to be written to a different folder, you can pass in the folder path to the `log_folder` argument.
 
+We're now ready to run the app. To do this we'll call the `fdmr::interactive_priors` function, passing in
+
+- `spatial_data = sp_data` - the spatial data covering areas across the city of Bristol
+- `measurement_data = covid19_data` - the COVID-19 data we want to model
+- `mesh = mesh` - the mesh we created above
+- `time_variable = week` - the name of the time variable in the data
+
 ```{r eval=FALSE}
 fdmr::interactive_priors(spatial_data = sp_data, measurement_data = covid19_data, mesh=mesh, time_variable="week")
 ```
@@ -86,6 +105,6 @@ fdmr::interactive_priors(spatial_data = sp_data, measurement_data = covid19_data
 
 The outputs of the model run by the app can be viewed in the `fdmr/logs` folder in your home directory. There you'll find three log files
 
-1. `priors_exploration_applog_timestamp.txt` - holds general logging information and errors
-2. `priors_exploration_parameters_timestamp.json` - holds the priors used in each model run
-3. `priors_exploration_modelout_timestamp.rds` - holds the model output for each run
+- `priors_exploration_applog_timestamp.txt` - holds general logging information and errors
+- `priors_exploration_parameters_timestamp.json` - holds the priors used in each model run
+- `priors_exploration_modelout_timestamp.rds` - holds the model output for each run