Skip to contents

This article is based on the blog post Integrating data from weathercan written for ROpenSci March 6th 2018.

In that article I demonstrated how we can incorporate data from weathercan into spatial visualizations. In this article I’d like to take that even further and show you how you can create interactive maps which highlight spatial variability in weather data.

Here, we’ll take a look at annual temperatures throughout different Eco Regions in Manitoba, Canada.

Setup

Using extra CSS styles

The map that we create here may look different for you unless you include the tweaks to the CSS styles I have made. If you’re using RMarkdown, you can supply these as a custom .css file, or inline with the <style> and </style> tags (like below).

<style>
div.leaflet-popup-content-wrapper {
  width: 700px;
  height: 100%;
}

div.leaflet-popup-tip-container {
  opacity: 0;
}

div.leaflet-popup-content {
  width: 90% !important;
}

img {
  max-width: 90% !important;
  min-width: 90% !important;
  border: none;
}

/* Fix the NA mis-alignment in the legend */
div.info.legend.leaflet-control br {
  clear: both;
}
</style>

Download Manitoban Eco Regions shapefile

download.file("http://mli2.gov.mb.ca/environment/shp_zip_files/env_ecological_areas_py_shp.zip",
              destfile = "./mapping_shp/ecological_shp.zip")
unzip("./mapping_shp/ecological_shp.zip")
file.remove("./mapping_shp/ecological_shp.zip")

Download Manitoba weather data

We’ll select all currently operating stations (end >= 2018) and download the daily weather data for 2017:

mb <- filter(stations(), prov == "MB", interval == "day", end >= 2018)
w <- weather_dl(mb$station_id, start = "2017-01-01", end = "2017-12-31", interval = "day")

Calculating summaries

In this section, we’ll summarize our data and create the means of showcasing it in our map. We’ll calculate some basic summaries and we’ll style some popups to contain this information (including the figures we’ll create later).

We’ll do a station-specific summaries/pop-ups for when users click on a station marker, and region-specific ones for when they click on the region polygon.

Station summaries

mb_stations <- w %>%
  group_by(station_id) %>%
  mutate(n = n(),
         n_missing = sum(is.na(mean_temp)),
         mean_temp = mean(mean_temp, na.rm = TRUE)) %>%
  filter((n - n_missing) > 0) %>% # Only keep stations with some temperature data
  select(station_name, station_id, lat, lon, mean_temp, n, n_missing) %>%
  distinct() %>%
  mutate(station_name = tools::toTitleCase(tolower(station_name)),
         info = paste0("<h3>Station: ", station_name, " (", station_id, ")</h2>",
                       "<hr>",
                       "<div>",
                       "<strong>Mean Temperature: </strong>", round(mean_temp, 1), "C<br>",
                       "<strong>No. days with data:  </strong>", n-n_missing, "<br>",
                       "<strong>No. days total:  </strong>", n, "</div>",
                       "<img src = './mapping_svg/", station_id, ".svg'>"),
         pretty_name = map(station_name, 
                           ~HTML(paste0("<strong>Station: </strong>", .x)))) %>%
  ungroup() %>%
  st_as_sf(coords = c("lon", "lat"), crs = "+proj=longlat")

Eco Region summaries

Before we summarize this data, we’ll filter the Eco Regions to just Manitoba and will join this to the stations data (after transforming the stations data to the same CRS), so we can figure out which stations belong to which regions.

mb_ecoregions <- st_read("./mapping_shp/env_ecological_areas.shp") %>%
  filter(MANITOBA == "yes") %>%
  # Get the larger-scale regions and combine
  group_by(ECOREGION, REGION_NAM, REGION_NOM) %>%
  summarize()
## Reading layer `env_ecological_areas' from data source 
##   `/tmp/RtmpU3VaDc/weathercan-source/vignettes/articles/mapping_shp/env_ecological_areas.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 279 features and 15 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -621621.3 ymin: 5386414 xmax: 1868189 ymax: 7166472
## Projected CRS: NAD83 / UTM zone 14N
## `summarise()` has grouped output by 'ECOREGION', 'REGION_NAM'. You can override
## using the `.groups` argument.
mb_stations <- st_transform(mb_stations, crs = st_crs(mb_ecoregions))
  
mb_ecoregions <- st_join(mb_ecoregions, mb_stations) %>%
  group_by(ECOREGION, REGION_NAM, REGION_NOM) %>%
  summarize(n_stations = length(unique(station_id[!is.na(station_id)])),
            mean_temp = mean(mean_temp, na.rm = TRUE),
            mean_temp = replace(mean_temp, is.nan(mean_temp), NA)) %>%
  mutate(info = paste0("<h3>Region: ", REGION_NAM, "/", REGION_NOM, " (", ECOREGION, ")</h2>",
                       "<hr>",
                       "<div>",
                       "<strong>Mean Temperature: </strong>", round(mean_temp, 1), 
                       if_else(!is.na(mean_temp), "C<br>", "<br>"),
                       "<strong>No. stations:  </strong>", n_stations, "</div>"),
         info = if_else(n_stations > 0, 
                        paste0(info, "<img src = './mapping_svg/", ECOREGION, ".svg'>"),
                        info),
         pretty_name = map(REGION_NAM, 
                           ~HTML(paste0("<strong>Region: </strong>", .x)))) %>%
  ungroup()
## `summarise()` has grouped output by 'ECOREGION', 'REGION_NAM'. You can override
## using the `.groups` argument.

Figures

We’ll set up some functions to create and save figures for our map

plot_station_fig <- function(d, station_id) {
  g <- ggplot(d, aes(x = date, y = mean_temp)) +
    theme_bw() +
    geom_line(na.rm = TRUE) +
    labs(x = "Date", y = "Mean Daily Temperature (C)")
  ggsave(paste0("./mapping_svg/", station_id, ".svg"), plot = g,
         width = 6, height = 3, dpi = 100)
}

plot_region_fig <- function(d, region) {
  g <- ggplot(d, aes(x = date, y = mean_temp, 
                group = station_name, colour = station_name)) +
    theme_bw() +
    geom_line(na.rm = TRUE) +
    scale_colour_viridis_d(end = 0.8) +
    labs(x = "Date", y = "Mean Daily Temperature (C)", colour = "Stations")
  ggsave(paste0("./mapping_svg/", region, ".svg"), plot = g,
         width = 6, height = 3, dpi = 100)
}

Now we’ll apply these functions to our data. Note that this figs object isn’t important, it’s just used as a way to loop through the data and save the figures as svg files.

figs <- st_join(mb_stations, mb_ecoregions) %>%
  st_set_geometry(NULL) %>%
  left_join(select(w, station_id, date, mean_temp), by = "station_id") %>%
  nest(-ECOREGION) %>%
  mutate(fig_region = map2(data, ECOREGION, ~plot_region_fig(.x, .y))) %>%
  unnest(data) %>%
  nest(-ECOREGION, -station_id) %>%
  mutate(fig_station = map2(data, station_id, ~plot_station_fig(.x, .y)))
## Warning: All elements of `...` must be named.
## Did you want `data = c(-ECOREGION, -station_id)`?
## Warning: All elements of `...` must be named.
## Did you want `data = -ECOREGION`?

Mapping

Finally, we’re ready to create our map!

We’ll start by transforming our data into WGS84 for leaflet, then we’ll create a palette and get some icons…

# Required for Leaflet
mb_ecoregions <- st_transform(mb_ecoregions, crs = 4326)
mb_stations <- st_transform(mb_stations, crs = 4326)

# Setup Palette for polygons
pal_eco <- colorNumeric(palette = "viridis",
                          domain = mb_ecoregions$mean_temp)

# Get icons for stations (red if above average, blue if below)
station_icons <- awesomeIcons(iconColor = "black", library = "ion",
                              markerColor = ifelse(mb_stations$mean_temp > 
                                                     mean(mb_stations$mean_temp, 
                                                          na.rm = TRUE), 
                                                   "red", "blue"))

Now for the real magic!

leaflet(width = "750px", height = "85vh") %>% 
  addTiles() %>%
  addPolygons(data = mb_ecoregions,
              color = "#444444", weight = 1, opacity = 1, fillOpacity = 0.5,
              fillColor = ~pal_eco(mean_temp),
              label = ~pretty_name, popup = ~info,
              popupOptions = popupOptions(keepInView = TRUE),
              highlightOptions = highlightOptions(bringToFront = TRUE, 
                                                  fillOpacity = 1)) %>%
  addAwesomeMarkers(data = mb_stations, group = "Stations",
                    icon = station_icons,
                    label = ~pretty_name, popup = ~info,
                    popupOptions = popupOptions(keepInView = TRUE)) %>%
  addLegend("bottomright", pal = pal_eco,
            values = mb_ecoregions$mean_temp,
            title = "Mean region temperature",
            labFormat = labelFormat(suffix = " C")) %>%
  addLegend("bottomright", 
            title = "Station temperature", 
            colors = c("#d63e2a", "#37a7da"), opacity = 1, 
            labels = c("Above the average", "Below the average")) %>%
  addLayersControl(
    overlayGroups = "Stations",
    options = layersControlOptions(collapsed = FALSE))

I think this is an interesting way of looking at the Eco Regions in Manitoba.

First we see a clear (and expected) pattern of decreasing temperatures with latitude (South to North). Further, we also see a bit of a South West to North East pattern.

Looking at the Eco Regions like this gives a clear idea of some of the parameters that make these Eco Regions distinct. For example, the patch of lands in the lower left corner of the province are the Mid-Boreal Uplands and the Boreal Transition. See how much cooler they are (on average) than the rest of southern Manitoba.

These southern Boreal regions represent a distinct ecology in southern Manitoba. If you zoom in on the map, you’ll see that they also hold two parks: Duck Mountain Provincial Park and Riding Mountain National Park.