diff --git a/CHANGELOG.md b/CHANGELOG.md
index f18fd8cd..7ebd9558 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -5,7 +5,25 @@ All notable changes to `fdmr` will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [Unreleased](https://github.com/openghg/openghg/compare/0.1.1...HEAD)
+## [Unreleased](https://github.com/openghg/openghg/compare/0.2.0...HEAD)
+
+## [0.2.0] - 2023-12-19
+
+### Fixed
+
+- Added a function to ensure correct polygon display across the dateline - [PR #285](https://github.com/4DModeller/fdmr/pull/285)
+
+### Added
+
+- Added a new tutorial on loading data from different sources - [PR #266](https://github.com/4DModeller/fdmr/pull/266/)
+- Added new geophysical processes tutorial thanks to Alexander Minakov (4minakov) - [PR #257](https://github.com/4DModeller/fdmr/pull/257)
+- Added mouse pointer coordinates header and standard measurement tool - [PR #260](https://github.com/4DModeller/fdmr/pull/260)
+- Added new help explainer to the `Help` tab of the `model_viewer` Shiny app - [PR #295](https://github.com/4DModeller/fdmr/pull/295)
+
+### Changed
+
+- Updated `plot_map` to allow use of both `leaflet` and `mapview` packages - [#291](https://github.com/4DModeller/fdmr/pull/291)
+- Moved from using `leaflet` to using (`mapview`)[https://r-spatial.github.io/mapview/index.html] for plotting the mesh and spatial data in the `mesh_builder` Shiny app. This enables use of UTM coordinates - [PR #288](https://github.com/4DModeller/fdmr/pull/288)
 
 ## [0.1.1] - 2023-11-01
 
@@ -25,6 +43,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Allowed reversal of colour palette and added raster plot legend in map plotter - [PR #229](https://github.com/4DModeller/fdmr/pull/229)
 - Added markers option to mesh plotter - [PR #230](https://github.com/4DModeller/fdmr/pull/230)
 - Updated the interfaces of the Shiny apps to the use [bslib](https://rstudio.github.io/bslib/index.html) theming - [PR #236](https://github.com/4DModeller/fdmr/pull/236)
+- Code tab added to `fdmr::model_viewer` so the user can easily reproduce plots - [PR #237](https://github.com/4DModeller/fdmr/pull/237)
 
 ## [0.1.0] - 2023-10-17
 
diff --git a/DESCRIPTION b/DESCRIPTION
index e2dcf469..c654b6b7 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: fdmr
 Title: 4D Modeller project
-Version: 0.1.1
+Version: 0.2.0
 Authors@R: 
     c(person("Gareth", "Jones", , "g.m.jones@bristol.ac.uk", role = c("aut", "cre"),
             comment = ""),
@@ -55,7 +55,12 @@ Imports:
     fmesher,
     purrr,
     shinyjs,
-    bslib
+    bslib,
+    bsicons,
+    future,
+    leafem,
+    sfheaders,
+    mapview
 Suggests: 
     bookdown,
     knitr,
diff --git a/NAMESPACE b/NAMESPACE
index fb7e01c5..dcf9a937 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,6 +1,8 @@
 # Generated by roxygen2: do not edit by hand
 
+export(antimeridian_wrapping)
 export(clear_caches)
+export(convert_from_lon_360)
 export(create_prediction_field)
 export(get_tutorial_datapath)
 export(latlong_to_utm)
diff --git a/R/plot_mapping.R b/R/plot_mapping.R
index c5ebb625..4ea0a1ba 100644
--- a/R/plot_mapping.R
+++ b/R/plot_mapping.R
@@ -1,4 +1,10 @@
-#' Create a simple Leaflet map from data
+#' Create a simple map from data using either the leaflet or mapview packages
+#'
+#' NOTE that the mapview backend is only intended for quick viewing of data,
+#' most of the customisation arguments are not available.
+#'
+#' The leaflet backend will work for most use cases and is recommended.
+#' For plotting of maps with UTM coordinates, the mapview backend is recommended.
 #'
 #' @param polygon_data Polygon data
 #' @param raster_data Raster datas
@@ -12,8 +18,11 @@
 #' @param polygon_line_colour Polygon surrounding line colour
 #' @param polygon_line_weight Polygon surrounding line weight
 #' @param reverse Reverse the colour palette if TRUE
+#' @param wrapping Split polygons along the antimeridian (-180/180 boundary) if TRUE
+#' @param backend Backend package to use for plotting, either "leaflet" or "mapview"
+#'
+#' @return leaflet::leaflet or mapview::mapview
 #'
-#' @return leaflet::leaflet
 #' @export
 plot_map <- function(polygon_data = NULL,
                      raster_data = NULL,
@@ -26,20 +35,84 @@ plot_map <- function(polygon_data = NULL,
                      polygon_line_colour = "grey",
                      polygon_line_weight = 1,
                      polygon_fill_opacity = 0.6,
-                     reverse = FALSE) {
+                     reverse = FALSE,
+                     wrapping = FALSE,
+                     backend = "leaflet") {
   if (is.null(polygon_data) && is.null(raster_data)) {
     stop("Polygon or raster data must be given.")
   }
-  library(leaflet)
+
+  if (backend == "leaflet") {
+    plot_map_leaflet(
+      polygon_data = polygon_data,
+      raster_data = raster_data,
+      domain = domain,
+      markers = markers,
+      palette = palette,
+      legend_title = legend_title,
+      add_scale_bar = add_scale_bar,
+      polygon_fill_colour = polygon_fill_colour,
+      polygon_line_colour = polygon_line_colour,
+      polygon_line_weight = polygon_line_weight,
+      polygon_fill_opacity = polygon_fill_opacity,
+      reverse = reverse,
+      wrapping = wrapping
+    )
+  } else if (backend == "mapview") {
+    plot_map_mapview(
+      polygon_data = polygon_data,
+      raster_data = raster_data
+    )
+  } else {
+    stop("Invalid backend given, must be either 'leaflet' or 'mapview'.")
+  }
+}
+
+
+#' Create a simple Leaflet map from data
+#'
+#' @param polygon_data Polygon data
+#' @param raster_data Raster data
+#' @param domain Domain data to be passed to leaflet::colorNumeric and leaflet::addLegend
+#' @param markers Markers to display on map. A named list with latitude, longitude and label names must be given.
+#' @param palette Palette to be used for colours, defaults to viridis
+#' @param legend_title Title for legend
+#' @param add_scale_bar Add scale bar if TRUE
+#' @param polygon_fill_opacity Leaflet polygon fill opacity, float from 0 to 1.0, passed to fillOpacity of leaflet::addPolygons
+#' @param polygon_fill_colour Polygon fill colour
+#' @param polygon_line_colour Polygon surrounding line colour
+#' @param polygon_line_weight Polygon surrounding line weight
+#' @param reverse Reverse the colour palette if TRUE
+#' @param wrapping Split polygons along the antimeridian (-180/180 boundary) if TRUE
+#'
+#' @return leaflet::leaflet
+#' @keywords internal
+plot_map_leaflet <- function(polygon_data = NULL,
+                             raster_data = NULL,
+                             domain = NULL,
+                             markers = NULL,
+                             palette = "viridis",
+                             legend_title = NULL,
+                             add_scale_bar = FALSE,
+                             polygon_fill_colour = "#E4572E",
+                             polygon_line_colour = "grey",
+                             polygon_line_weight = 1,
+                             polygon_fill_opacity = 0.6,
+                             reverse = FALSE,
+                             wrapping = FALSE) {
   m <- leaflet::leaflet()
   m <- leaflet::addTiles(m)
-  m <- leaflet::addProviderTiles(m, leaflet::providers$Esri.WorldImagery, group = "Satellite")
+  m <- leaflet::addProviderTiles(m, leaflet::providers$Openstreetmap, group = "Satellite")
+  m <- leafem::addMouseCoordinates(m, native.crs = TRUE)
 
   # Store a vector of layers we add to the map,
   # used later to create the layers control object
   layers <- c()
 
   if (!is.null(polygon_data)) {
+    if (isTRUE(wrapping)) {
+      polygon_data <- fdmr::antimeridian_wrapping(polygon_data, crs = "+proj=longlat +datum=WGS84", unique_inst = TRUE, to_sp = FALSE)
+    }
     if (!is.null(domain)) {
       colours <- leaflet::colorNumeric(palette = palette, domain = domain, reverse = reverse)
       polygon_fill_colour <- ~ colours(domain)
@@ -63,7 +136,7 @@ plot_map <- function(polygon_data = NULL,
   }
 
   if (!is.null(raster_data)) {
-    colours <- leaflet::colorNumeric(palette = palette, domain = raster::values(raster_data), na.color=rgb(0,0,0,0), reverse = reverse)
+    colours <- leaflet::colorNumeric(palette = palette, domain = raster::values(raster_data), na.color = rgb(0, 0, 0, 0), reverse = reverse)
     m <- leaflet::addRasterImage(m,
       x = raster_data,
       opacity = 0.75,
@@ -72,11 +145,11 @@ plot_map <- function(polygon_data = NULL,
       colors = colours,
     )
     m <- leaflet::addLegend(m,
-                            pal = colours,
-                            values = raster::values(raster_data),
-                            opacity = 0.75,
-                            title = legend_title,
-                            na.label = ""
+      pal = colours,
+      values = raster::values(raster_data),
+      opacity = 0.75,
+      title = legend_title,
+      na.label = ""
     )
     layers <- append(layers, "Raster")
   }
@@ -97,5 +170,28 @@ plot_map <- function(polygon_data = NULL,
     m <- leaflet::addScaleBar(m, position = "bottomleft")
   }
 
+  m <- leaflet::addMeasure(m, position = "bottomleft", primaryLengthUnit = "kilometers", primaryAreaUnit = "sqmeters")
+
+  return(m)
+}
+
+
+#' A simple map plotter using mapview. This is only intended for very quick viewing of data.
+#'
+#' @param polygon_data Spatial data to plot
+#' @param raster_data Raster data to plot
+#'
+#' @return mapview::mapview
+#' @keywords internal
+plot_map_mapview <- function(polygon_data = NULL, raster_data = NULL) {
+  map_types <- c("OpenStreetMap", "Esri.WorldImagery", "OpenTopoMap")
+
+  m <- mapview::mapview(map.types = map_types)
+  if (!is.null(polygon_data)) {
+    m <- m + mapview::mapview(polygon_data)
+  }
+  if (!is.null(raster_data)) {
+    m <- m + mapview::mapview(raster_data)
+  }
   return(m)
 }
diff --git a/R/plot_mesh.R b/R/plot_mesh.R
index f83144ca..617ad9a4 100644
--- a/R/plot_mesh.R
+++ b/R/plot_mesh.R
@@ -24,7 +24,9 @@ plot_mesh <- function(mesh, spatial_data = NULL, longitude_column = "LONG", lati
     )
   }
 
-  spatial_mesh <- fdmr::mesh_to_spatial(mesh = mesh, crs = expected_crs)
+  spatial_mesh_original <- fdmr::mesh_to_spatial(mesh = mesh, crs = expected_crs)
+
+  spatial_mesh <- fdmr::antimeridian_wrapping(spatial_mesh_original, crs = expected_crs, unique_inst = FALSE, to_sp = FALSE)
 
   plot_polygons <- FALSE
   plot_points <- FALSE
@@ -60,6 +62,8 @@ plot_mesh <- function(mesh, spatial_data = NULL, longitude_column = "LONG", lati
   m <- leaflet::leaflet()
   m <- leaflet::addTiles(m, group = "OSM")
   m <- leaflet::addPolygons(m, data = spatial_mesh, weight = 0.5, fillOpacity = 0.2, fillColor = "#5252ca", group = "Mesh")
+  m <- leaflet::addMeasure(m, position = "bottomleft", primaryLengthUnit = 'kilometers', primaryAreaUnit = 'sqmeters')
+  m <- leafem::addMouseCoordinates(m, native.crs = TRUE)
 
   if (plot_polygons) {
     m <- leaflet::addPolygons(m, data = spatial_data, fillColor = "#d66363", color = "green", weight = 1, group = "Spatial")
diff --git a/R/plot_priors.R b/R/plot_priors.R
index 6258df31..7c8524a0 100644
--- a/R/plot_priors.R
+++ b/R/plot_priors.R
@@ -4,10 +4,9 @@
 #' @param to_plot Type of data to plot, "Range for f" etc
 #'
 #' @return ggplot2::ggplot
-#' @keywords internal
 plot_line_comparison <- function(data, to_plot, title) {
-  ar1_data <- purrr::map(data, function(x) as.data.frame(x$pars[[to_plot]]))
-  single_df <- dplyr::bind_rows(ar1_data, .id = "Run")
+  parsed_data <- purrr::map(data, function(x) as.data.frame(x$pars[[to_plot]]))
+  single_df <- dplyr::bind_rows(parsed_data, .id = "Run")
   if (nrow(single_df) == 0) {
     return("No pars data.")
   }
@@ -18,32 +17,11 @@ plot_line_comparison <- function(data, to_plot, title) {
     ggplot2::theme(text = ggplot2::element_text(size = 16))
 }
 
-
-#' Plot AR(1)
-#'
-#' @param data Parsed model output
-#' @param to_plot Type of data to plot, "Range for f" etc
-#'
-#' @return ggplot2::ggplot
-#' @keywords internal
-plot_ar1 <- function(data) {
-  ar1_data <- purrr::map(data, function(x) as.data.frame(x$pars$`GroupRho for f`))
-  single_df <- dplyr::bind_rows(ar1_data, .id = "Run")
-  if (nrow(single_df) == 0) {
-    return("No pars data.")
-  }
-
-  ggplot2::ggplot(single_df, ggplot2::aes(x = x, y = y, color = Run)) +
-    ggplot2::geom_line() +
-    ggplot2::theme(text = ggplot2::element_text(size = 16))
-}
-
 #' Create boxplots from priors run data
 #'
 #' @param data
 #'
 #' @return graphics::boxplot
-#' @keywords internal
 plot_priors_boxplot <- function(data) {
   # TODO - I'm sure this can be done in a nicer functional way
   fitted_mean_post <- purrr::map(data, function(x) x$fitted_mean_post)
@@ -60,7 +38,6 @@ plot_priors_boxplot <- function(data) {
 #' @param measurement_data Measurement data
 #'
 #' @return ggplot2::ggplot
-#' @keywords internal
 plot_priors_density <- function(data, measurement_data) {
   # Can this be done in a cleaner way? Just create a dataframe from the lists?
   fitted_values <- unlist(purrr::map(data, function(x) x$fitted_mean_post))
@@ -82,7 +59,6 @@ plot_priors_density <- function(data, measurement_data) {
 #' @param data
 #'
 #' @return ggplot2::ggplot
-#' @keywords internal
 plot_dic <- function(data) {
   infocri <- base::cbind.data.frame(
     priors = unlist(purrr::map(seq(1, length(data)), function(x) paste("Run", x))),
diff --git a/R/shiny_meshbuilder.R b/R/shiny_meshbuilder.R
index d72711d4..8a22be0d 100644
--- a/R/shiny_meshbuilder.R
+++ b/R/shiny_meshbuilder.R
@@ -1,13 +1,13 @@
 #' Mesh building shiny app
 #'
-#' @param spatial_data a data.frame or tibble containing spatial data
-#' @param data Observations data, for use with the check_mesh functionality
+#' @param spatial_data Spatial data
+#' @param obs_data Measurement data
 #' @param crs CRS as a proj4string
 #' @param offset Specifies the size of the inner and outer extensions around data locations, passed to fmesher::fm_mesh_2d_inla
 #' @param max_edge The largest allowed triangle edge length. One or two values, passed to fmesher::fm_mesh_2d_inla
 #' @param cutoff The minimum allowed distance between points, passed to fmesher::fm_mesh_2d_inla
-#' @param latitude_column Name of the latitude column in the spatial data
-#' @param longitude_column Name of the longitude column in the spatial data
+#' @param y_coord Name of the latitude column in the spatial data
+#' @param x_coord Name of the longitude column in the spatial data
 #'
 #' @importFrom magrittr %>%
 #'
@@ -17,12 +17,12 @@ meshbuilder_shiny <- function(
     spatial_data,
     obs_data = NULL,
     crs = NULL,
-    max_edge = NULL,
     offset = NULL,
+    max_edge = NULL,
     cutoff = NULL,
     plot_poly = FALSE,
-    latitude_column = "LAT",
-    longitude_column = "LONG") {
+    y_coord = "LAT",
+    x_coord = "LONG") {
   if (!is.data.frame(spatial_data) && !is(spatial_data, "SpatialPolygonsDataFrame") && !is(spatial_data, "SpatialPointsDataFrame")) {
     stop("spatial_data must be a data.frame or tibble containing columns with latitude and longitude data.")
   }
@@ -47,81 +47,95 @@ meshbuilder_shiny <- function(
     )
   }
 
-  got_lat_long <- all(c(longitude_column, latitude_column) %in% names(spatial_data))
+  # If the user passes in any of these then we enable the sliders
+  enable_inputs <- (!is.null(max_edge) || !is.null(offset) || !is.null(cutoff))
+
+  got_lat_long <- all(c(x_coord, y_coord) %in% names(spatial_data))
   if (!got_lat_long) {
     stop("Cannot read latitude and longitude data from spatial data. Please ensure given names are correct.")
   }
 
-  # The number of nodes we count as being a big mesh
-  n_nodes_big_mesh <- 10000
-
-  default_max_edge <- c(0.1, 0.3)
-  default_offset <- c(0.2, 0.7)
-  default_cutoff <- 0.2
+  default_max_edge_min <- 0.01
+  default_max_edge_max <- 0.3
+  default_offset_min <- 0.02
+  default_offset_max <- 0.2
+  default_cutoff <- 0.02
   # TODO - these defaults need changing?
-  if (is.null(max_edge)) max_edge <- default_max_edge
-  if (is.null(offset)) offset <- default_offset
-  if (is.null(cutoff)) cutoff <- default_cutoff
+  if (!is.null(max_edge)) {
+    default_max_edge_min <- max_edge[1]
+    default_max_edge_max <- max_edge[2]
+  }
+
+  if (!is.null(offset)) {
+    default_offset_min <- offset[1]
+    default_offset_max <- offset[2]
+  }
+
+  if (!is.null(cutoff)) default_cutoff <- cutoff
 
   # Make sure we have our own internal correctly formatted version of the data
-  coords_only <- spatial_data[, c(longitude_column, latitude_column)]
+  coords_only <- spatial_data[, c(x_coord, y_coord)]
   names(coords_only) <- c("LONG", "LAT")
 
-  plot_polygons <- FALSE
-  plot_points <- FALSE
-  # We may not use this
-  spatial_points <- NULL
-  # Do some checks to see what kind of data we have
-  # If we have a SpatialPolygonsDataFrame, we can plot the polygons
-  # otherwise if we just have SpatialPoints we can plot the points
-  # otherwise we don't plot anything
-  if (class(spatial_data) == "SpatialPolygonsDataFrame") {
-    plot_polygons <- TRUE
-  } else {
-    plot_points <- TRUE
-
-    spatial_points <- sp::SpatialPointsDataFrame(
-      coords = coords_only,
-      data = spatial_data,
-      proj4string = sp::CRS(crs)
-    )
-  }
 
   busy_spinner <- get_busy_spinner()
 
-  ui <- bslib::page_fluid(
-    theme = bslib::bs_theme(bootswatch = "cosmo"),
+  ui <- shiny::fluidPage(
     busy_spinner,
-    shinybusy::add_loading_state(
-      "#map",
-      timeout = 600,
-      text = "Calculating mesh...",
-      svgColor = "steelblue"
-    ),
+    shinyjs::useShinyjs(),
     shiny::headerPanel(title = "Creating a mesh"),
     shiny::sidebarLayout(
       shiny::sidebarPanel(
-        shiny::sliderInput(
-          inputId = "max_edge",
-          label = "Max edge:",
-          min = 0.02, value = c(0.1, 0.3), max = 10
+        shiny::checkboxInput(inputId = "enable_inputs", label = "Enable customisation", value = enable_inputs),
+        shiny::h4("Max edge"),
+        shiny::fluidRow(
+          shiny::column(
+            6,
+            shiny::numericInput(
+              inputId = "max_edge_min",
+              label = "Min:",
+              value = default_max_edge_min
+            )
+          ),
+          shiny::column(
+            6,
+            shiny::numericInput(
+              inputId = "max_edge_max",
+              label = "Max:",
+              value = default_max_edge_max
+            )
+          )
         ),
         shiny::p("Max permitted edge length for a triangle"),
-        shiny::sliderInput(
-          inputId = "offset",
-          label = "Offset:",
-          min = 0.02, value = c(0.2, 0.7), max = 10
+        shiny::h4("Offset"),
+        shiny::fluidRow(
+          shiny::column(
+            6,
+            shiny::numericInput(
+              inputId = "offset_min",
+              label = "Min:",
+              value = default_offset_min
+            ),
+          ),
+          shiny::column(
+            6,
+            shiny::numericInput(
+              inputId = "offset_max",
+              label = "Max:",
+              value = default_offset_max
+            )
+          )
         ),
         shiny::p("Specifies the size of the inner and outer extensions around data locations."),
-        shiny::sliderInput(
+        shiny::h4("Cutoff"),
+        shiny::numericInput(
           inputId = "cutoff",
-          label = "Cutoff:",
-          min = 0.005, value = 0.2, max = 0.9
+          label = NULL,
+          value = default_cutoff
         ),
         shiny::p("Minimum allowed distance between data points."),
         shiny::actionButton("plot_mesh", label = "Plot mesh"),
         shiny::actionButton("reset_mesh", label = "Reset"),
-        # shiny::actionButton("check_button", "Check mesh"),
       ),
       shiny::mainPanel(
         shiny::tabsetPanel(
@@ -129,7 +143,9 @@ meshbuilder_shiny <- function(
           shiny::tabPanel(
             "Plot",
             class = "p-3 border",
-            shiny::div(id = "map_div", leaflet::leafletOutput("map", height = "80vh"))
+            shiny::div(leaflet::leafletOutput("map", height = "80vh")),
+            shiny::br(),
+            shiny::textOutput(outputId = "mesh_crs")
           ),
           shiny::tabPanel("Code", class = "p-3 border", shiny::verbatimTextOutput("mesh_code")),
           shiny::tabPanel(
@@ -151,92 +167,99 @@ meshbuilder_shiny <- function(
 
   # Define server logic required to draw a histogram
   server <- function(input, output, session) {
+    shiny::observeEvent(input$enable_inputs, {
+      if (input$enable_inputs) {
+        shinyjs::enable("max_edge_min")
+        shinyjs::enable("max_edge_max")
+        shinyjs::enable("offset_min")
+        shinyjs::enable("offset_max")
+        shinyjs::enable("cutoff")
+      } else {
+        shinyjs::disable("max_edge_min")
+        shinyjs::disable("max_edge_max")
+        shinyjs::disable("offset_min")
+        shinyjs::disable("offset_max")
+        shinyjs::disable("cutoff")
+      }
+    })
+
     shiny::observeEvent(input$reset_mesh, {
-      shiny::updateSliderInput(session, inputId = "max_edge", value = default_max_edge)
-      shiny::updateSliderInput(session, inputId = "offset", value = default_offset)
-      shiny::updateSliderInput(session, inputId = "cutoff", value = default_cutoff)
+      shiny::updateNumericInput(session, inputId = "max_edge_min", value = default_max_edge_min)
+      shiny::updateNumericInput(session, inputId = "max_edge_max", value = default_max_edge_max)
+      shiny::updateNumericInput(session, inputId = "offset_min", value = default_offset_min)
+      shiny::updateNumericInput(session, inputId = "offset_max", value = default_offset_max)
+      shiny::updateNumericInput(session, inputId = "cutoff", value = default_cutoff)
+    })
+
+    output$mesh_crs <- shiny::renderText({
+      paste("Mesh CRS: ", crs)
     })
 
     mesh <- shiny::eventReactive(input$plot_mesh, ignoreNULL = FALSE, {
+      if (input$enable_inputs) {
+        max_edge <- c(input$max_edge_min, input$max_edge_max)
+        offset <- c(input$offset_min, input$offset_max)
+        cutoff <- input$cutoff
+      } else {
+        max_edge <- NULL
+        offset <- NULL
+        cutoff <- NULL
+      }
+
       fmesher::fm_mesh_2d_inla(
         loc = coords_only,
-        max.edge = input$max_edge,
-        cutoff = input$cutoff,
-        offset = input$offset,
+        max.edge = max_edge,
+        cutoff = cutoff,
+        offset = offset,
         crs = crs,
       )
     })
 
-    # large_mesh <- shiny::reactive({
-    #   mesh()$n > n_nodes_big_mesh
-    # })
-
     mesh_spatial <- shiny::reactive(
       suppressMessages(
         suppressWarnings(
-          fdmr::mesh_to_spatial(mesh = mesh(), crs = crs)
+          fdmr::antimeridian_wrapping(fdmr::mesh_to_spatial(mesh = mesh(), crs = crs), crs = crs, unique_inst = FALSE, to_sp = FALSE)
         )
       )
     )
 
-    output$map <- leaflet::renderLeaflet({
-      m <- leaflet::leaflet()
-      m <- leaflet::addTiles(m, group = "OSM")
-      m <- leaflet::addPolygons(m, data = mesh_spatial(), weight = 0.5, fillOpacity = 0.2, fillColor = "#5252ca", group = "Mesh")
-      if (plot_polygons) {
-        m <- leaflet::addPolygons(m, data = spatial_data, fillColor = "#d66363", color = "green", weight = 1, group = "Spatial")
-      } else if (plot_points) {
-        m <- leaflet::addCircles(m, data = spatial_points, group = "Spatial", fillColor = "#b9220b", color = "#b9220b")
+    spatial <- shiny::reactive({
+      if (is.data.frame(spatial_data)) {
+        sf::st_as_sf(
+          spatial_data,
+          coords = c(x_coord, y_coord),
+          crs = crs
+        )
+      } else {
+        spatial_data
       }
-
-      m <- leaflet::addLayersControl(m,
-        position = "topright",
-        baseGroups = c("OSM"),
-        overlayGroups = c("Mesh", "Spatial"),
-        options = leaflet::layersControlOptions(collapsed = FALSE)
-      )
     })
 
-    output$mesh_code <- shiny::reactive(
-      paste0(
-        "location_data <- spatial_data[, c('", longitude_column, "', '", latitude_column, "')]\n\n",
-        "names(location_data) <- c('LONG', 'LAT')\n\n",
-        "mesh <- fmesher::fm_mesh_2d_inla(loc = location_data,
-          max.edge = c(", paste0(input$max_edge, collapse = ", "), "),
-          cutoff = ", input$cutoff, ",
-          offset=c(", paste0(input$offset, collapse = ", "), "))\n"
-      )
-    )
-
-    # shiny::observe({
-    #   if (large_mesh() && !modal_shown()) {
-    #     shiny::showModal(shiny::modalDialog(
-    #       "Mesh is large, plotting may be slow.",
-    #       title = "Mesh warning",
-    #       easyClose = TRUE,
-    #       footer = NULL
-    #     ))
-    #   }
-    #   modal_shown(TRUE)
-    # })
-
+    output$map <- leaflet::renderLeaflet({
+      map_tiles <- c("OpenStreetMap", "Esri.WorldImagery", "OpenTopoMap")
+      m <- mapview::mapview(mesh_spatial(), layer.name = "Mesh", col.regions = "#548C2F", map.types = map_tiles) + mapview::mapview(spatial(), layer.name = "Spatial")
+      m@map
+    })
 
-    shiny::observeEvent(input$check_button, {
-      if (is.null(obs_data) || is.null(mesh())) {
-        errors <- "No observation data. Cannot check mesh."
+    output$mesh_code <- shiny::reactive({
+      if (input$enable_inputs) {
+        max_edge_str <- paste0("max.edge = c(", input$max_edge_min, ",", input$max_edge_max, "),")
+        cutoff_str <- paste0("cutoff = ", input$cutoff, ",")
+        offset_str <- paste0("offset = c(", input$offset_min, ",", input$max_edge_max, ")")
       } else {
-        errors <- fdmr::mesh_checker(mesh = mesh(), observations = obs_data)
-        if (!length(errors)) {
-          errors <- "No errors found."
-        }
+        max_edge_str <- "max.edge = NULL,"
+        cutoff_str <- "cutoff = NULL,"
+        offset_str <- "offset = NULL"
       }
 
-      shiny::showModal(shiny::modalDialog(
-        stringr::str_flatten(errors, collapse = "\n"),
-        title = "Mesh check",
-        easyClose = TRUE,
-        footer = NULL
-      ))
+      paste0(
+        "location_data <- spatial_data[, c('", x_coord, "', '", y_coord, "')],\n",
+        "names(location_data) <- c('LONG', 'LAT')\n",
+        "mesh <- fmesher::fm_mesh_2d_inla(loc = location_data,\n\t",
+        max_edge_str, "\n\t",
+        cutoff_str, "\n\t",
+        offset_str, ")\n"
+      )
     })
   }
 
@@ -253,12 +276,12 @@ meshbuilder_shiny <- function(
 #' @param offset Specifies the size of the inner and outer extensions around data locations, passed to fmesher::fm_mesh_2d_inla
 #' @param max_edge The largest allowed triangle edge length. One or two values, passed to fmesher::fm_mesh_2d_inla
 #' @param cutoff The minimum allowed distance between points, passed to fmesher::fm_mesh_2d_inla
-#' @param latitude_column Name of the latitude column in the spatial data
-#' @param longitude_column Name of the longitude column in the spatial data
+#' @param y_coord Name of the latitude column in the spatial data
+#' @param x_coord Name of the longitude column in the spatial data
 #'
 #' @return shiny::app
 #' @export
-mesh_builder <- function(spatial_data, obs_data = NULL, crs = NULL, max_edge = NULL, offset = NULL, cutoff = NULL, latitude_column = "LAT", longitude_column = "LONG") {
+mesh_builder <- function(spatial_data, obs_data = NULL, crs = NULL, max_edge = NULL, offset = NULL, cutoff = NULL, y_coord = "LAT", x_coord = "LONG") {
   shiny::runApp(meshbuilder_shiny(
     spatial_data = spatial_data,
     obs_data = obs_data,
@@ -266,7 +289,7 @@ mesh_builder <- function(spatial_data, obs_data = NULL, crs = NULL, max_edge = N
     max_edge = max_edge,
     offset = offset,
     cutoff = cutoff,
-    latitude_column = latitude_column,
-    longitude_column = longitude_column
+    y_coord = y_coord,
+    x_coord = x_coord
   ))
 }
diff --git a/R/shiny_modelbuilder.R b/R/shiny_modelbuilder.R
index 34b4b6da..fc395d31 100644
--- a/R/shiny_modelbuilder.R
+++ b/R/shiny_modelbuilder.R
@@ -604,7 +604,9 @@ model_builder_shiny <- function(spatial_data,
       leaflet::leaflet() %>%
         leaflet::addTiles(group = "OSM") %>%
         leaflet::addRasterImage(map_raster(), colors = map_colours(), opacity = 0.9, group = "Raster") %>%
-        leaflet::addLegend(position = "topright", pal = map_colours(), values = legend_values())
+        leaflet::addLegend(position = "topright", pal = map_colours(), values = legend_values()) %>%
+        leaflet::addMeasure(position = "bottomleft", primaryLengthUnit = 'kilometers', primaryAreaUnit = 'sqmeters') %>%
+        leafem::addMouseCoordinates(native.crs = TRUE)
     })
 
     model_plot <- shiny::eventReactive(input$plot_type, ignoreNULL = FALSE, {
diff --git a/R/shiny_modelviewer.R b/R/shiny_modelviewer.R
index 23580c30..25be16e7 100644
--- a/R/shiny_modelviewer.R
+++ b/R/shiny_modelviewer.R
@@ -45,12 +45,14 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data, data_distri
       type = "tabs",
       shiny::tabPanel(
         "Plots",
+        class = "p-3 border",
         shiny::h2("Plot output"),
         shiny::selectInput(inputId = "plot_type", label = "Plot type:", choices = plot_choices, selected = plot_choices[1]),
         shiny::plotOutput(outputId = "plot_model_out")
       ),
       shiny::tabPanel(
         "Map",
+        class = "p-3 border",
         shiny::fluidRow(
           shiny::column(
             4,
@@ -80,9 +82,29 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data, data_distri
         ),
         leaflet::leafletOutput(outputId = "map_out")
       ),
+      shiny::tabPanel(
+        "Code",
+        class = "p-3 border",
+        shiny::verbatimTextOutput(outputId = "code_out")
+      ),
       shiny::tabPanel(
         "Help",
-        shiny::h3("Help"),
+        class = "p-3 border",
+        shiny::h4("Range"),
+        shiny::p("The Range plot provides a plot of the posterior distribution of the spatial range parameter ρ, the distance at which the correlation between two points is approximately 0."),
+        shiny::h4("Stdev"),
+        shiny::p("The Stdev plot provides a plot of the posterior distribution of the marginal standard deviation σ of the spatial field."),
+        shiny::h4("AR(1)"),
+        shiny::p("The AR(1) plot provides a plot of the posterior distribution of the temporal autocorrelation α. A value close to 1 indicates strong positive temporal dependence, a value of 0 indicates independence across time, and a value close to -1 indicates strong negative temporal dependence."),
+        shiny::h4("Boxplot"),
+        shiny::p("The Boxplot displays the predicted (or fitted) values across the observed locations and time points, showing how the predictions are distributed."),
+        shiny::h4("Density"),
+        shiny::p("The Density plot displays the density curve of the predicted (or fitted) values across the observed locations and time points."),
+        shiny::h4("DIC"),
+        shiny::p("The DIC plot displays the Deviance Information Criterion (DIC) value for the model. DIC is a model selection criterion, with lower values indicating better model fit."),
+        shiny::h4("Map"),
+        shiny::p("The plot of predicted mean fields maps the predicted values across a number of grid locations that cover the study region. It helps to identify areas with high predictions and those with low predictions. The predicted values at observed locations and time points can be obtained by using the syntax model_output$summary.fitted.values$mean[1:nrow(observed_data)]. The predicted values at the grid locations are calculated using the method described in the priors tutorial "),
+        shiny::p("The plot of random effect fields maps the random effect values (the f() values) across a number of grid locations that cover the study region. The random effect values at the mesh nodes and time points can be obtained by using the syntax model_output$summary.random$f$mean.")
       )
     )
   )
@@ -203,6 +225,63 @@ model_viewer_shiny <- function(model_output, mesh, measurement_data, data_distri
       }
     })
 
+    output$code_out <- shiny::reactive({
+      code_str <- ""
+      parsed_model_str <- "parsed_model_out <- fdmr::parse_model_output(model_output = model_output,
+                                                      measurement_data = measurement_data)"
+
+      if (input$plot_type %in% c("Range", "Stdev", "AR(1)")) {
+        if (input$plot_type == "Range") {
+          to_plot <- "Range for f"
+          title <- "Range"
+        } else if (input$plot_type == "Stdev") {
+          to_plot <- "Stdev for f"
+          title <- "Marginal standard deviation"
+        } else if (input$plot_type == "AR(1)") {
+          to_plot <- "GroupRho for f"
+          title <- "AR(1)"
+        }
+
+        code_str <- paste0('parsed_data <- purrr::map(parsed_model_out, function(x) as.data.frame(x$pars[["', to_plot, '"]]))
+        single_df <- dplyr::bind_rows(parsed_data, .id = "Run")
+
+        ggplot2::ggplot(single_df, ggplot2::aes(x = x, y = y, color = Run)) +
+        ggplot2::geom_line() +
+        ggplot2::ggtitle("', title, '") +
+        ggplot2::theme(text = ggplot2::element_text(size = 16))')
+      } else if (input$plot_type == "Boxplot") {
+        code_str <- 'fitted_mean_post <- purrr::map(parsed_model_out, function(x) x$fitted_mean_post)
+        names(fitted_mean_post) <- purrr::map(seq(1, length(parsed_model_out)), function(x) paste("Run", x))
+        post_rate <- cbind.data.frame(fitted_mean_post)
+        graphics::boxplot(post_rate, xlab = "Prior scenario", ylab = "Fitted values")'
+      } else if (input$plot_type == "Density") {
+        code_str <- 'fitted_values <- unlist(purrr::map(parsed_model_out, function(x) x$fitted_mean_post))
+          run_strings <- unlist(purrr::map(seq(1, length(parsed_model_out)), function(x) paste("Run", x)))
+
+          post_rate <- base::cbind.data.frame(
+            "Prior scenario" = rep(run_strings, each = nrow(measurement_data)),
+            "Fitted values" = fitted_values
+          )
+
+          ggplot2::ggplot(post_rate, ggplot2::aes(x = `Fitted values`, color = `Prior scenario`)) +
+            ggplot2::geom_density() +
+            ggplot2::theme(text = ggplot2::element_text(size = 16))'
+      } else if (input$plot_type == "DIC") {
+        code_str <- 'infocri <- base::cbind.data.frame(
+        priors = unlist(purrr::map(seq(1, length(parsed_model_out)), function(x) paste("Run", x))),
+        DIC = unlist(purrr::map(parsed_model_out, function(x) x$dic))
+      )
+
+      infocri$priors <- base::as.factor(infocri$priors)
+
+      ggplot2::ggplot(infocri, ggplot2::aes(x = priors, y = DIC)) +
+        ggplot2::geom_point() +
+        ggplot2::theme(text = ggplot2::element_text(size = 16))'
+      }
+
+      paste0(parsed_model_str, "\n\n", code_str)
+    })
+
     output$plot_model_out <- shiny::renderPlot({
       model_plot()
     })
diff --git a/R/util_coords.R b/R/util_coords.R
index 07c7a649..def37fa6 100644
--- a/R/util_coords.R
+++ b/R/util_coords.R
@@ -40,3 +40,112 @@ has_coords <- function(spatial_data) {
     }
   )
 }
+
+#' Convert longitudes from 0 to 360 degrees to -180 to 180 degrees
+#'
+#' @param sf_data An sf object; does not accept SpatialPolygon* objects
+#' @param crs CRS as a proj4string or EPSG code
+#' @param add_data Select if data associated with the object are carried forward by the transformed version, defaults to FALSE
+#' @param longitude_column Name of longitude, defaults to LONG
+#'
+#' @return polygons with converted coordinates
+#' @export
+convert_from_lon_360 <- function(sf_data, crs = 4326, add_data = TRUE, longitude_column = "LONG") {
+  # Get polygon coordinates
+  coords <- sf::st_coordinates(sf_data)
+  
+  # Check if POLYGON contains holes, remove if so
+  if ((base::length(base::unique(coords[, "L1"])) > 1) & (base::ncol(coords) > 3)) {
+    warning("Converter does not define the interior holes. Please convert your data yourself if these features are crucial")
+    sf_data <- sfheaders::sf_remove_holes(obj = sf_data)
+    coords <- sf::st_coordinates(sf_data)
+  }
+  
+  coords_conv <- coords
+  
+  # Convert data coordinates if selected
+  if (base::isTRUE(add_data) & (!base::is.null(longitude_column))) {
+    if (base::is.null(sf_data[["LONG"]]) & (longitude_column == "LONG")) {
+      warning("`LONG` not found. Please change the default `longitude_column` name to the longitude name in your dataset")
+    } else {
+      sf_data[[longitude_column]] <- base::ifelse(sf_data[[longitude_column]] > 180, 
+                                                  sf_data[[longitude_column]] - 360, 
+                                                  sf_data[[longitude_column]])
+    }
+  }
+  
+  # Convert polygon longitudes
+  coords_conv[, "X"] <- base::ifelse(coords[, "X"] > 180, coords[, "X"] - 360, coords[, "X"])
+  
+  # Define the feature column and convert to POLYGON
+  if (base::ncol(coords_conv) == 4) {
+    poly <- "L2"
+    conv <- sfheaders::sf_polygon(obj = coords_conv, x = 'X', y = 'Y', polygon_id = poly)
+    conv <- sf::st_sf(conv, crs = crs)
+  } else if (base::ncol(coords_conv) == 5) {
+    poly <- "L3"
+    multipoly <- "L2"
+    conv <- sfheaders::sf_multipolygon(obj = coords_conv, x = 'X', y = 'Y', multipolygon_id = multipoly, polygon_id = poly)
+    conv <- sf::st_sf(conv, crs = crs)
+  } else {
+    stop("Object type not supported (must be POLYGON or MULTIPOLYGON)")
+  }
+  
+  # Add the associated data if requested
+  if (base::isTRUE(add_data)) {
+    r_conv <- sf_data
+    r_conv$geometry <- conv$geometry
+  } else {
+    r_conv <- conv
+  }
+  
+  # Return a list of POLYGON objects with shifted longitudes
+  return(r_conv)
+}
+
+#' Split polygons into two if crossing the dateline (antimeridian)
+#' 
+#' @param polygon_data SpatialPolygon or SpatialPolygonDataFrame object
+#' @param crs CRS as a proj4string or EPSG code, default EPSG:4326 (datum=WGS84)
+#' @param unique_inst Option to remove duplicate polygons (with the same coordinates), default TRUE
+#' @param to_sp Option to convert POLYGON object into SpatialPolygon*, default TRUE
+
+#' @return SpatialPolygon or SpatialPolygonDataFrame split along the antimeridian
+#' @export
+antimeridian_wrapping <- function(polygon_data, crs = 4326, unique_inst = TRUE, to_sp = TRUE) {
+  # Convert to sf object if needed
+  if (!is(polygon_data, 'sf')) {
+    if (base::isTRUE(unique_inst)) {
+      all_coords <- base::lapply(polygon_data@polygons, function(x) {
+        base::lapply(x@Polygons, function(x) {
+          sp::coordinates(x)
+        })
+      })
+      polygon_data <- polygon_data[!base::duplicated(all_coords),]
+    }
+    sf_data <- sf::st_as_sf(polygon_data)
+    if (!isTRUE(sf::st_is_longlat(sf_data))) {
+      stop("Projected coordinates not allowed, please choose a different projection to avoid the antimeridian problem")
+    }
+  } else {
+    sf_data <- polygon_data
+  }
+  
+  # Check if longitudes go from -180 to 180 rather than 0 to 360
+  if (sf::st_bbox(sf_data)[['xmin']] > 0) {
+    polygon_data_conv <- convert_from_lon_360(sf_data, crs = crs, add_data = FALSE)
+    warning("Polygon coordinates [0;360] have been converted to [-180;180]")
+  }
+  
+  # Apply the function
+  sf_split <- sf::st_wrap_dateline(sf_data, options = c("WRAPDATELINE=YES", "DATELINEOFFSET=90"))
+  
+  # Revert sf objects back to sp
+  if (base::isTRUE(to_sp)) {
+    polygon_data_conv <- sf::as_Spatial(sf_split)
+  } else {
+    polygon_data_conv <- sf_split
+  }
+  
+  return(polygon_data_conv)
+}
diff --git a/README.md b/README.md
index 3b1176fd..7a25c950 100644
--- a/README.md
+++ b/README.md
@@ -20,7 +20,7 @@
 
 ### - evaluate the fully trained model output
 
-![](https://github.com/4DModeller/fdmr/assets/8915182/fe791f74-c9b4-4db0-a52e-ffa018b12b41)
+![](https://github.com/4DModeller/fdmr/assets/8915182/f2b68ccc-cfba-4393-b2f4-6cecdad332c0)
 
 ## Quickstart
 
@@ -58,4 +58,14 @@ sudo apt-get install libharfbuzz-dev libfribidi-dev libfreetype6-dev \
                         libpng-dev libtiff5-dev libjpeg-dev libudunits2-dev libgdal-dev
 ```
 
-Note that on other Linux distributions the names of these packages may differ.
\ No newline at end of file
+Note that on other Linux distributions the names of these packages may differ.
+
+## Want to contribute?
+
+You can contribute to 4DModeller in a variety of ways including: responding to issues, introducing new features such as new tutorials or core functionality, or helping to plan a future 4DModeller hackathon. See below how to do each:
+
+1. *Issues:* Please [checkout our issues page](https://github.com/4DModeller/fdmr/issues). If you see something you can solve then [fork the repo, make the changes, then make a pull request](https://stackoverflow.com/questions/14587045/how-to-merge-branch-of-forked-repo-into-master-branch-of-original-repo#14587354). If you have an issue with 4DModeller, please open an issue instead.
+2. *New Features:* new features can be handled in two ways. First, you can suggest new features [using the GitHub issue tracker](https://github.com/4DModeller/fdmr/issues). Second, you can contribute new features by [forking the repo, creating the new tutorial or core functionality, then making a pull request](https://github.com/4DModeller/fdmr/issues).
+3. *Hackathon Planning:* If you would like to help organize a 4DModeller hackathon either by helping organize a core hackathon or by organizing one yourself at your institution, then please reach out to one of the 4DModeller developers.
+
+If you make regular contributions through issues and new features then we would be happy to include you in the core group as a developer of 4DModeller.
diff --git a/_pkgdown.yml b/_pkgdown.yml
index cfebb6a8..4abdf5c0 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -86,10 +86,12 @@ reference:
   desc: Functions to help convert coordinates and units
   contents:
     - latlong_to_utm
-    - to_dates
+    - antimeridian_wrapping
+    - convert_from_lon_360
 - title: Utility functions
   desc:  These help us handle paths, retrieve data etc
   contents:
+    - to_dates
     - retrieve_tutorial_data
     - get_tutorial_datapath
     - load_tutorial_data
diff --git a/man/antimeridian_wrapping.Rd b/man/antimeridian_wrapping.Rd
new file mode 100644
index 00000000..2eaf85e4
--- /dev/null
+++ b/man/antimeridian_wrapping.Rd
@@ -0,0 +1,28 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/util_coords.R
+\name{antimeridian_wrapping}
+\alias{antimeridian_wrapping}
+\title{Split polygons into two if crossing the dateline (antimeridian)}
+\usage{
+antimeridian_wrapping(
+  polygon_data,
+  crs = 4326,
+  unique_inst = TRUE,
+  to_sp = TRUE
+)
+}
+\arguments{
+\item{polygon_data}{SpatialPolygon or SpatialPolygonDataFrame object}
+
+\item{crs}{CRS as a proj4string or EPSG code, default EPSG:4326 (datum=WGS84)}
+
+\item{unique_inst}{Option to remove duplicate polygons (with the same coordinates), default TRUE}
+
+\item{to_sp}{Option to convert POLYGON object into SpatialPolygon*, default TRUE}
+}
+\value{
+SpatialPolygon or SpatialPolygonDataFrame split along the antimeridian
+}
+\description{
+Split polygons into two if crossing the dateline (antimeridian)
+}
diff --git a/man/convert_from_lon_360.Rd b/man/convert_from_lon_360.Rd
new file mode 100644
index 00000000..0bb32900
--- /dev/null
+++ b/man/convert_from_lon_360.Rd
@@ -0,0 +1,28 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/util_coords.R
+\name{convert_from_lon_360}
+\alias{convert_from_lon_360}
+\title{Convert longitudes from 0 to 360 degrees to -180 to 180 degrees}
+\usage{
+convert_from_lon_360(
+  sf_data,
+  crs = 4326,
+  add_data = TRUE,
+  longitude_column = "LONG"
+)
+}
+\arguments{
+\item{sf_data}{An sf object; does not accept SpatialPolygon* objects}
+
+\item{crs}{CRS as a proj4string or EPSG code}
+
+\item{add_data}{Select if data associated with the object are carried forward by the transformed version, defaults to FALSE}
+
+\item{longitude_column}{Name of longitude, defaults to LONG}
+}
+\value{
+polygons with converted coordinates
+}
+\description{
+Convert longitudes from 0 to 360 degrees to -180 to 180 degrees
+}
diff --git a/man/mesh_builder.Rd b/man/mesh_builder.Rd
index 30ca5dde..b8b00eee 100644
--- a/man/mesh_builder.Rd
+++ b/man/mesh_builder.Rd
@@ -11,8 +11,8 @@ mesh_builder(
   max_edge = NULL,
   offset = NULL,
   cutoff = NULL,
-  latitude_column = "LAT",
-  longitude_column = "LONG"
+  y_coord = "LAT",
+  x_coord = "LONG"
 )
 }
 \arguments{
@@ -28,9 +28,9 @@ mesh_builder(
 
 \item{cutoff}{The minimum allowed distance between points, passed to fmesher::fm_mesh_2d_inla}
 
-\item{latitude_column}{Name of the latitude column in the spatial data}
+\item{y_coord}{Name of the latitude column in the spatial data}
 
-\item{longitude_column}{Name of the longitude column in the spatial data}
+\item{x_coord}{Name of the longitude column in the spatial data}
 }
 \value{
 shiny::app
diff --git a/man/meshbuilder_shiny.Rd b/man/meshbuilder_shiny.Rd
index 673598a1..27656e22 100644
--- a/man/meshbuilder_shiny.Rd
+++ b/man/meshbuilder_shiny.Rd
@@ -8,30 +8,30 @@ meshbuilder_shiny(
   spatial_data,
   obs_data = NULL,
   crs = NULL,
-  max_edge = NULL,
   offset = NULL,
+  max_edge = NULL,
   cutoff = NULL,
   plot_poly = FALSE,
-  latitude_column = "LAT",
-  longitude_column = "LONG"
+  y_coord = "LAT",
+  x_coord = "LONG"
 )
 }
 \arguments{
-\item{spatial_data}{a data.frame or tibble containing spatial data}
+\item{spatial_data}{Spatial data}
 
-\item{crs}{CRS as a proj4string}
+\item{obs_data}{Measurement data}
 
-\item{max_edge}{The largest allowed triangle edge length. One or two values, passed to fmesher::fm_mesh_2d_inla}
+\item{crs}{CRS as a proj4string}
 
 \item{offset}{Specifies the size of the inner and outer extensions around data locations, passed to fmesher::fm_mesh_2d_inla}
 
-\item{cutoff}{The minimum allowed distance between points, passed to fmesher::fm_mesh_2d_inla}
+\item{max_edge}{The largest allowed triangle edge length. One or two values, passed to fmesher::fm_mesh_2d_inla}
 
-\item{latitude_column}{Name of the latitude column in the spatial data}
+\item{cutoff}{The minimum allowed distance between points, passed to fmesher::fm_mesh_2d_inla}
 
-\item{longitude_column}{Name of the longitude column in the spatial data}
+\item{y_coord}{Name of the latitude column in the spatial data}
 
-\item{data}{Observations data, for use with the check_mesh functionality}
+\item{x_coord}{Name of the longitude column in the spatial data}
 }
 \value{
 shiny::app
diff --git a/man/plot_ar1.Rd b/man/plot_ar1.Rd
deleted file mode 100644
index f5b26b62..00000000
--- a/man/plot_ar1.Rd
+++ /dev/null
@@ -1,20 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/plot_priors.R
-\name{plot_ar1}
-\alias{plot_ar1}
-\title{Plot AR(1)}
-\usage{
-plot_ar1(data)
-}
-\arguments{
-\item{data}{Parsed model output}
-
-\item{to_plot}{Type of data to plot, "Range for f" etc}
-}
-\value{
-ggplot2::ggplot
-}
-\description{
-Plot AR(1)
-}
-\keyword{internal}
diff --git a/man/plot_dic.Rd b/man/plot_dic.Rd
index 69fff308..b92fa7d9 100644
--- a/man/plot_dic.Rd
+++ b/man/plot_dic.Rd
@@ -15,4 +15,3 @@ ggplot2::ggplot
 \description{
 Plot Deviance Information Criterion (DIC) values
 }
-\keyword{internal}
diff --git a/man/plot_line_comparison.Rd b/man/plot_line_comparison.Rd
index f953e712..b185b7d4 100644
--- a/man/plot_line_comparison.Rd
+++ b/man/plot_line_comparison.Rd
@@ -17,4 +17,3 @@ ggplot2::ggplot
 \description{
 Plot line comparison for stdev etc
 }
-\keyword{internal}
diff --git a/man/plot_map.Rd b/man/plot_map.Rd
index 6daa45b9..75c2f3d9 100644
--- a/man/plot_map.Rd
+++ b/man/plot_map.Rd
@@ -2,7 +2,7 @@
 % Please edit documentation in R/plot_mapping.R
 \name{plot_map}
 \alias{plot_map}
-\title{Create a simple Leaflet map from data}
+\title{Create a simple map from data using either the leaflet or mapview packages}
 \usage{
 plot_map(
   polygon_data = NULL,
@@ -15,7 +15,10 @@ plot_map(
   polygon_fill_colour = "#E4572E",
   polygon_line_colour = "grey",
   polygon_line_weight = 1,
-  polygon_fill_opacity = 0.6
+  polygon_fill_opacity = 0.6,
+  reverse = FALSE,
+  wrapping = FALSE,
+  backend = "leaflet"
 )
 }
 \arguments{
@@ -40,10 +43,21 @@ plot_map(
 \item{polygon_line_weight}{Polygon surrounding line weight}
 
 \item{polygon_fill_opacity}{Leaflet polygon fill opacity, float from 0 to 1.0, passed to fillOpacity of leaflet::addPolygons}
+
+\item{reverse}{Reverse the colour palette if TRUE}
+
+\item{wrapping}{Split polygons along the antimeridian (-180/180 boundary) if TRUE}
+
+\item{backend}{Backend package to use for plotting, either "leaflet" or "mapview"}
 }
 \value{
-leaflet::leaflet
+leaflet::leaflet or mapview::mapview
 }
 \description{
-Create a simple Leaflet map from data
+NOTE that the mapview backend is only intended for quick viewing of data,
+most of the customisation arguments are not available.
+}
+\details{
+The leaflet backend will work for most use cases and is recommended.
+For plotting of maps with UTM coordinates, the mapview backend is recommended.
 }
diff --git a/man/plot_map_leaflet.Rd b/man/plot_map_leaflet.Rd
new file mode 100644
index 00000000..b6e583c5
--- /dev/null
+++ b/man/plot_map_leaflet.Rd
@@ -0,0 +1,56 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_mapping.R
+\name{plot_map_leaflet}
+\alias{plot_map_leaflet}
+\title{Create a simple Leaflet map from data}
+\usage{
+plot_map_leaflet(
+  polygon_data = NULL,
+  raster_data = NULL,
+  domain = NULL,
+  markers = NULL,
+  palette = "viridis",
+  legend_title = NULL,
+  add_scale_bar = FALSE,
+  polygon_fill_colour = "#E4572E",
+  polygon_line_colour = "grey",
+  polygon_line_weight = 1,
+  polygon_fill_opacity = 0.6,
+  reverse = FALSE,
+  wrapping = FALSE
+)
+}
+\arguments{
+\item{polygon_data}{Polygon data}
+
+\item{raster_data}{Raster data}
+
+\item{domain}{Domain data to be passed to leaflet::colorNumeric and leaflet::addLegend}
+
+\item{markers}{Markers to display on map. A named list with latitude, longitude and label names must be given.}
+
+\item{palette}{Palette to be used for colours, defaults to viridis}
+
+\item{legend_title}{Title for legend}
+
+\item{add_scale_bar}{Add scale bar if TRUE}
+
+\item{polygon_fill_colour}{Polygon fill colour}
+
+\item{polygon_line_colour}{Polygon surrounding line colour}
+
+\item{polygon_line_weight}{Polygon surrounding line weight}
+
+\item{polygon_fill_opacity}{Leaflet polygon fill opacity, float from 0 to 1.0, passed to fillOpacity of leaflet::addPolygons}
+
+\item{reverse}{Reverse the colour palette if TRUE}
+
+\item{wrapping}{Split polygons along the antimeridian (-180/180 boundary) if TRUE}
+}
+\value{
+leaflet::leaflet
+}
+\description{
+Create a simple Leaflet map from data
+}
+\keyword{internal}
diff --git a/man/plot_map_mapview.Rd b/man/plot_map_mapview.Rd
new file mode 100644
index 00000000..85e024e4
--- /dev/null
+++ b/man/plot_map_mapview.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_mapping.R
+\name{plot_map_mapview}
+\alias{plot_map_mapview}
+\title{A simple map plotter using mapview. This is only intended for very quick viewing of data.}
+\usage{
+plot_map_mapview(polygon_data = NULL, raster_data = NULL)
+}
+\arguments{
+\item{polygon_data}{Spatial data to plot}
+
+\item{raster_data}{Raster data to plot}
+}
+\value{
+mapview::mapview
+}
+\description{
+A simple map plotter using mapview. This is only intended for very quick viewing of data.
+}
+\keyword{internal}
diff --git a/man/plot_mesh.Rd b/man/plot_mesh.Rd
index 928db33e..9656973d 100644
--- a/man/plot_mesh.Rd
+++ b/man/plot_mesh.Rd
@@ -8,7 +8,8 @@ plot_mesh(
   mesh,
   spatial_data = NULL,
   longitude_column = "LONG",
-  latitude_column = "LAT"
+  latitude_column = "LAT",
+  markers = NULL
 )
 }
 \arguments{
@@ -20,6 +21,9 @@ that can be converted to a data.frame with longitude and latitude columns}
 \item{longitude_column}{Longitude column in spatial_data}
 
 \item{latitude_column}{Latitude column in spatial_data name}
+
+\item{markers}{Markers to display on top of mesh. A named list with latitude, longitude and label names must be given.
+Expects longitude name to be longitude, latitude name to be latitude, label name to be label.}
 }
 \value{
 leaflet::leaflet
diff --git a/man/plot_priors_boxplot.Rd b/man/plot_priors_boxplot.Rd
index 2ca7e7cb..e7ee38dc 100644
--- a/man/plot_priors_boxplot.Rd
+++ b/man/plot_priors_boxplot.Rd
@@ -15,4 +15,3 @@ graphics::boxplot
 \description{
 Create boxplots from priors run data
 }
-\keyword{internal}
diff --git a/man/plot_priors_density.Rd b/man/plot_priors_density.Rd
index ef702e46..c840ab5c 100644
--- a/man/plot_priors_density.Rd
+++ b/man/plot_priors_density.Rd
@@ -17,4 +17,3 @@ ggplot2::ggplot
 \description{
 Plot density function
 }
-\keyword{internal}
diff --git a/vignettes/covid.Rmd b/vignettes/covid.Rmd
index c6be9b43..23123147 100644
--- a/vignettes/covid.Rmd
+++ b/vignettes/covid.Rmd
@@ -4,7 +4,7 @@ output:
   bookdown::html_document2:
     base_format: rmarkdown::html_vignette
     fig_caption: yes
-bibliography: "covid.bib"
+bibliography: "references.bib"
 link-citations: yes
 vignette: >
   %\VignetteIndexEntry{COVID-19}
diff --git a/vignettes/covid_ST_mcmc.Rmd b/vignettes/covid_ST_mcmc.Rmd
index 2acdb33d..003c1fe9 100644
--- a/vignettes/covid_ST_mcmc.Rmd
+++ b/vignettes/covid_ST_mcmc.Rmd
@@ -4,7 +4,7 @@ output:
   bookdown::html_document2:
     base_format: rmarkdown::html_vignette
     fig_caption: yes
-bibliography: "covid.bib"
+bibliography: "references.bib"
 link-citations: yes
 vignette: >
   %\VignetteIndexEntry{COVID-19 MCMC (Space-time)}
diff --git a/vignettes/covid_mcmc.Rmd b/vignettes/covid_mcmc.Rmd
index 95d1e832..e86e16c5 100644
--- a/vignettes/covid_mcmc.Rmd
+++ b/vignettes/covid_mcmc.Rmd
@@ -4,7 +4,7 @@ output:
   bookdown::html_document2:
     base_format: rmarkdown::html_vignette
     fig_caption: yes
-bibliography: "covid.bib"
+bibliography: "references.bib"
 link-citations: yes
 vignette: >
   %\VignetteIndexEntry{COVID-19 MCMC}
@@ -218,7 +218,7 @@ sp_data@data$est.rate.mcmc <- mcmc_fitted_prev
 domain <- sp_data@data$est.rate.mcmc
 legend_values <- sp_data@data$est.rate.mcmc
 
-fdmr::plot_map(
+fdmr::plot_map_leaflet(
   polygon_data = sp_data,
   domain = domain,
   palette = "Reds",
@@ -236,7 +236,7 @@ sp_data@data$est.rate.inlabym <- INLA_model$summary.fitted.values$mean
 domain <- sp_data@data$est.rate.inlabym
 legend_values <- sp_data@data$est.rate.inlabym
 
-fdmr::plot_map(
+fdmr::plot_map_leaflet(
   polygon_data = sp_data,
   domain = domain,
   palette = "Reds",
@@ -254,7 +254,7 @@ sp_data@data$est.rate.inlaspde <- inlabru_model$summary.fitted.values$mean[1:nro
 domain <- sp_data@data$est.rate.inlaspde
 legend_values <- sp_data@data$est.rate.inlaspde
 
-fdmr::plot_map(
+fdmr::plot_map_leaflet(
   polygon_data = sp_data,
   domain = domain,
   palette = "Reds",
diff --git a/vignettes/data/MT_Svalbard_Z.csv b/vignettes/data/MT_Svalbard_Z.csv
new file mode 100644
index 00000000..ca612830
--- /dev/null
+++ b/vignettes/data/MT_Svalbard_Z.csv
@@ -0,0 +1,37 @@
+T,Code,Lat,Lon,X,Y,Zxx,DZxx
+64,W01,79.439116,13.371895,344352.8496,8829185.331,0.52156-1.2687i,0.054863
+128,W01,79.439116,13.371895,344352.8496,8829185.331,0.97673-0.94615i,0.079632
+256,W01,79.439116,13.371895,344352.8496,8829185.331,1.275-0.54561i,0.10554
+64,W02,79.471702,13.403403,345465.1801,8832708.491,0.36581+0.027603i,0.012205
+128,W02,79.471702,13.403403,345465.1801,8832708.491,0.3344+0.13025i,0.019636
+256,W02,79.471702,13.403403,345465.1801,8832708.491,0.33298+0.071971i,0.016023
+64,W03,79.35887,14.091742,357911.5152,8818456.864,0.67767+0.53303i,0.01449
+128,W03,79.35887,14.091742,357911.5152,8818456.864,0.6028+0.52425i,0.022896
+256,W03,79.35887,14.091742,357911.5152,8818456.864,0.46475+0.40616i,0.029995
+64,W04,79.392861,14.0056,356602.2036,8822435.623,-0.44539-1.2946i,0.012704
+128,W04,79.392861,14.0056,356602.2036,8822435.623,-0.12597-0.99152i,0.016897
+256,W04,79.392861,14.0056,356602.2036,8822435.623,0.13967-0.52094i,0.021547
+64,W05,79.372686,13.864863,353455.5873,8820549.936,0.8007+0.63546i,0.018563
+128,W05,79.372686,13.864863,353455.5873,8820549.936,0.90669+0.78722i,0.025465
+256,W05,79.372686,13.864863,353455.5873,8820549.936,0.54578+0.61109i,0.022635
+64,W06,79.337495,13.893527,353563.6793,8816578.44,0.85439+0.74466i,0.039878
+128,W06,79.337495,13.893527,353563.6793,8816578.44,1.0516+1.1238i,0.054053
+256,W06,79.337495,13.893527,353563.6793,8816578.44,0.49092+0.83855i,0.042233
+64,W07,79.42051667,13.736,351485.8998,8826175.357,0.81914+0.80551i,0.023629
+128,W07,79.42051667,13.736,351485.8998,8826175.357,0.78128+0.68671i,0.030426
+256,W07,79.42051667,13.736,351485.8998,8826175.357,0.51883+0.59256i,0.030384
+64,W08,79.431243,13.699478,350892.5935,8827456.802,0.39831-0.1402i,0.01706
+128,W08,79.431243,13.699478,350892.5935,8827456.802,0.44428-0.09936i,0.022626
+256,W08,79.431243,13.699478,350892.5935,8827456.802,0.42242+0.15965i,0.028102
+64,W09,79.49686111,13.91783333,356218.9673,8834177.87,-0.12304-0.52841i,0.00536
+128,W09,79.49686111,13.91783333,356218.9673,8834177.87,-0.011371-0.37717i,0.0055408
+256,W09,79.49686111,13.91783333,356218.9673,8834177.87,0.040909-0.16656i,0.0060876
+64,W10,79.472456,13.944975,356438.2811,8831406.374,0.19744-0.74573i,0.014357
+128,W10,79.472456,13.944975,356438.2811,8831406.374,0.4019-0.45596i,0.021548
+256,W10,79.472456,13.944975,356438.2811,8831406.374,0.50063+0.068497i,0.027314
+64,W11,79.486643,13.235664,342293.8848,8834811.661,0.82246+2.7602i,0.016108
+128,W11,79.486643,13.235664,342293.8848,8834811.661,0.38066+1.8882i,0.02087
+256,W11,79.486643,13.235664,342293.8848,8834811.661,0.17339+1.0451i,0.027693
+64,W12,79.425669,13.378223,344285.4562,8827679.973,0.69293-0.38905i,0.032038
+128,W12,79.425669,13.378223,344285.4562,8827679.973,0.90861-0.030763i,0.0516
+256,W12,79.425669,13.378223,344285.4562,8827679.973,0.813+0.39348i,0.073448
diff --git a/vignettes/data_preprocessing.Rmd b/vignettes/data_preprocessing.Rmd
index 9403597c..3f373109 100644
--- a/vignettes/data_preprocessing.Rmd
+++ b/vignettes/data_preprocessing.Rmd
@@ -4,7 +4,7 @@ output:
   bookdown::html_document2:
     base_format: rmarkdown::html_vignette
     fig_caption: yes
-bibliography: "covid.bib"
+bibliography: "references.bib"
 link-citations: yes
 vignette: >
   %\VignetteIndexEntry{Data pre-processing}
diff --git a/vignettes/decision_tree.Rmd b/vignettes/decision_tree.Rmd
index 09a92854..7ea02596 100644
--- a/vignettes/decision_tree.Rmd
+++ b/vignettes/decision_tree.Rmd
@@ -107,11 +107,3 @@ fit <- bru(value ~ 1 + observed(x + y, model = "linear"), data = df)
 # Print the model summary
 print(summary(fit))
 ```
-
-### Comparing the output
-
-Finally, we can compare the output of all the models using the DIC. This is important as we would like to know the fittedness of each model even if we have some theoretical belief that a model should be specificied in a certain way.
-
-```{r error=TRUE}
-plot(DICs)
-```
diff --git a/vignettes/meshbuilder.Rmd b/vignettes/meshbuilder.Rmd
index 70bf7b92..7a029780 100644
--- a/vignettes/meshbuilder.Rmd
+++ b/vignettes/meshbuilder.Rmd
@@ -41,7 +41,7 @@ Let's call the `mesh_builder` function and pass in the data and tell it to use `
 > **_NOTE:_**  The mesh builder may take a short time (~ 20s) to build this mesh due to its size. 
 
 ```{r eval=FALSE}
-fdmr::mesh_builder(spatial_data = lakes_data, longitude_column = "centroid_lon", latitude_column = "centroid_lat")
+fdmr::mesh_builder(spatial_data = lakes_data, x_coord = "centroid_lon", y_coord = "centroid_lat")
 ``` 
 
 ## Loading data - `rds` file
@@ -81,29 +81,14 @@ offset <- c(initial_range / 4, initial_range)
 cutoff <- max_edge / 7
 ```
 
+> **_NOTE:_** By default we let `fmesher` select the defaults for these values itself. If you encounter long mesh build times try using the `meshbuilder` defaults of NULL values for `max_edge` etc. 
+
 Now we're ready to start the app.
 
 ```{r eval=FALSE}
 fdmr::mesh_builder(spatial_data = sp_data, max_edge = max_edge_fin, offset = offset, cutoff = cutoff)
 ```
 
-## Run checks on mesh
-
-We provide a simple function to check meshes created using the mesh builder tool. From with the user interface click "Check mesh" to run a number of tests on the mesh.
-Currently these are:
-
-1. Check that the number of mesh nodes isn't greater than the number of measurments
-2. Check that the number of triangles isn't greater than the number of measurements
-3. Check that there are no isolated triangles
-
-To use the mesh checking functionality you must pass your measuremnet data to the `fdmr::mesh_builder` function. Here we use the COVID-19 data. Create a mesh of your design and when you're finished click the "Check mesh" button. This passes the created mesh
-to the `fdmr::mesh_checker` function and returns a list containing any errors found with the mesh.
-
-```{r eval=FALSE}
-covid19_data <- fdmr::load_tutorial_data(dataset = "covid", filename = "covid19_data.rds")
-fdmr::mesh_builder(spatial_data = sp_data, obs_data = covid19_data)
-```
-
 ## Exporting your mesh
 
 To export your mesh click on the Code tab and copy and paste the code used to created the mesh.
diff --git a/vignettes/priors.Rmd b/vignettes/priors.Rmd
index 92c44c07..3f131a17 100644
--- a/vignettes/priors.Rmd
+++ b/vignettes/priors.Rmd
@@ -4,7 +4,7 @@ output:
   bookdown::html_document2:
     base_format: rmarkdown::html_vignette
     fig_caption: yes
-bibliography: "covid.bib"
+bibliography: "references.bib"
 link-citations: yes
 vignette: >
   %\VignetteIndexEntry{Priors exploration}
diff --git a/vignettes/covid.bib b/vignettes/references.bib
similarity index 100%
rename from vignettes/covid.bib
rename to vignettes/references.bib
diff --git a/vignettes/simulation_gaussian_data.Rmd b/vignettes/simulation_gaussian_data.Rmd
index a93efc6e..8e528241 100644
--- a/vignettes/simulation_gaussian_data.Rmd
+++ b/vignettes/simulation_gaussian_data.Rmd
@@ -225,11 +225,7 @@ fdmr::plot_mesh(mesh = mesh, spatial_data = simdf@coords)
 fdmr::mesh_builder(spatial_data = sp_data)
 ```
 
-The interactive map allows you to customise the initial parameters such as `max_edge`, `offset` and `cutoff` to change the shape and size of the mesh. We also provide a simple function to check if there are any errors with the meshes created using the mesh builder tool. To use the mesh checking functionality you must pass your measuremnet data to the `fdmr::mesh_builder` function. Specifically, create a mesh of your design using the meshbuilder tool. When you’re finished, click the “Check mesh” button on the user interface. This passes the created mesh to the `fdmr::mesh_checker` function and returns a list containing any errors found with the mesh.
-
-```{r meshchecking,eval=FALSE}
-fdmr::mesh_builder(spatial_data = sp_data, obs_data = simdf)
-```
+The interactive map allows you to customise the initial parameters such as `max_edge`, `offset` and `cutoff` to change the shape and size of the mesh.
 
 ## Model builder Shiny app
 `fdmr` provides a Priors Shiny app which allows you to interactively set and see the model fitting results of different priors. You can launch this app by passing the spatial and observation data to the `fdmr::model_builder` function.
diff --git a/vignettes/simulation_poisson_data.Rmd b/vignettes/simulation_poisson_data.Rmd
index b1656e63..c4f22e6b 100644
--- a/vignettes/simulation_poisson_data.Rmd
+++ b/vignettes/simulation_poisson_data.Rmd
@@ -227,11 +227,7 @@ fdmr::plot_mesh(mesh = mesh, spatial_data = simdf@coords)
 ```{r meshbuilder,eval=FALSE}
 fdmr::mesh_builder(spatial_data = sp_data)
 ```
-The interactive map allows you to customise the initial parameters such as `max_edge`, `offset` and `cutoff` to change the shape and size of the mesh. We also provide a simple function to check if there are any errors with the meshes created using the mesh builder tool. To use the mesh checking functionality you must pass your measuremnet data to the `fdmr::mesh_builder` function. Specifically, create a mesh of your design using the meshbuilder tool. When you’re finished, click the “Check mesh” button on the user interface. This passes the created mesh to the `fdmr::mesh_checker` function and returns a list containing any errors found with the mesh.
-
-```{r meshchecking,eval=FALSE}
-fdmr::mesh_builder(spatial_data = sp_data, obs_data = simdf)
-```
+The interactive map allows you to customise the initial parameters such as `max_edge`, `offset` and `cutoff` to change the shape and size of the mesh. 
 
 ## Model builder Shiny app
 `fdmr` provides a Priors Shiny app which allows you to interactively set and see the model fitting results of different priors. You can launch this app by passing the spatial and observation data to the `fdmr::model_builder` function.
diff --git a/vignettes/tutorial_geophysics.Rmd b/vignettes/tutorial_geophysics.Rmd
new file mode 100644
index 00000000..ccc31fea
--- /dev/null
+++ b/vignettes/tutorial_geophysics.Rmd
@@ -0,0 +1,106 @@
+---
+title: "4DModeller for geophysical signals"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{4DModeller for geophysical signals}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+## Introduction
+
+This is a tutorial to apply R-INLA to modeling geophysical data. The data represent time-series at some locations distributed over a line or over an area. There are two possible case studies one is seismic data with seismic stations located along a line. Each station record 4 components of acoustic signal X-Y-Z particle displacement on geophone and presure component on hydrophone sensor. The second case is magnteotelluric data with stations distrubuted over some area. The signal has 4 channels: 2 magnetic and 2 electric field channels. In the dataset only 2-3 stations measure at the same. The spatiotemporal evolution of the field is governed by Maxwell equations. The source of EM field is disturbance of ionosphere due to solar activity, electrical structure of the crust and noise component (cultural EM noise, wind, rain, local conductors). The goal is to describe the source signal components in space in time.  
+The spatiotemporal evolution of the field is governed by wave equation. The small-scale heterogenetities in earth crust produce multiply scattered wavefield, getting more expressed at a later times and called coda wave. The goal is to learn about correlations between signals at diffent stations and from this predict distribution of heterogeneities.     
+
+## Pre-processing and Import data
+
+Geophysical data can come in lots of different formats but here will will start with a csv.
+
+```{r}
+# Set the path to the CSV file
+data_path <- "data/MT_Svalbard_Z.csv"
+
+# Read the CSV file into a data frame
+d <- read.csv(data_path)
+
+# Display the first 30 rows of the data frame
+print(head(d, 30))
+
+# Create a scatter plot
+plot(d$X, d$Y, pch = 20, main = "Svalbard MT sites", xlab = "X", ylab = "Y")
+```
+
+Here we create the locations of the data for the mesh creation. This is the locations of the data, not the locations of the mesh nodes.
+```{r}
+locations <- d[, c("Lon", "Lat")]
+locations <- unique(locations)
+names(locations) <- c("LONG", "LAT")
+```
+
+Here we use the locations to help specify the best fit mesh to our data.
+
+## Meshing
+```{r}
+mesh <- fmesher::fm_mesh_2d(
+  loc.domain = locations,
+  max.edge = 0.05,
+  cutoff = 1e-3,
+  offset = 0.1
+)
+plot(mesh)
+points(locations, col = "red")
+
+fdmr::plot_mesh(mesh)
+```
+
+## Stochastic modeling with INLA
+
+Now we can create a model that is solved on the mesh.
+
+```{r}
+library(INLA)
+
+# Synthetic data generation
+set.seed(123) # For reproducibility
+n <- 100 # Number of time points
+stations <- 12 # Number of stations
+
+# Generate time index
+time_index <- 1:n
+
+# Generate spatial index (station IDs)
+space_index <- rep(1:stations, each = n)
+
+# Simulate some harmonic signals with noise for three stations
+harmonic_data <- data.frame(
+  station = factor(space_index),
+  time = rep(time_index, stations),
+  observation = sin(rep(time_index, stations) * 2 * pi / 50) +
+    rnorm(n * stations, sd = 0.5) +
+    rep(rnorm(stations, sd = 3), each = n) # Station-specific offset
+)
+
+# Define the model with harmonic terms for time and spatial correlation
+formula <- observation ~ f(station, model = "iid") +
+  f(time, model = "rw1", cyclic = TRUE)
+
+# Fit the model using INLA
+result <- inla(formula, family = "gaussian", data = harmonic_data)
+
+# Display the summary of the results
+summary(result)
+
+# Visualize the fitted values
+plot(harmonic_data$time, harmonic_data$observation, col = harmonic_data$station, pch = 19, cex = 0.5, xlab = "Time", ylab = "Observation", main = "MT signals at 12 stations")
+points(harmonic_data$time, result$summary.fitted.values$mean, pch = 4, cex = 0.7, col = "blue")
+legend("topright", legend = c("Observations", "Fitted"), col = c("black", "blue"), pch = c(19, 4))
+
+
+# Extract the hyperparameters of the spatial field
+spatial_hyperparams <- result$summary.hyperpar
+
+# Print the hyperparameters
+print(spatial_hyperparams)
+```
+
+## Observed time series will be added soon