In sustainablefarms/sflddata: Data Extraction Methods for the Linking Data Project

knitr::opts_chunk$set(echo = TRUE, root.dir = "..")
library(sflddata)

A document for storing / demonstrating methods for plotting RasterBricks with many layers.

invisible(lapply(c("raster", "maptools", "rgdal", "ncdf4", "lubridate", "sf"),
                 library, character.only = TRUE))
sws_sites <- readRDS("../private/data/clean/sws_sites.rds")
points <- sws_sites_2_spdf(sws_sites)
b <- brick_fmc(points, 2001)

Using ggplot2

library(ggplot2)
library(rasterVis)
library(viridis)
library(ggrepel)

Facet Grids using two different methods:

gplot(subset(b, 1:2)) +
  geom_tile(aes(fill = value)) +
  facet_grid(~ variable) + 
  scale_fill_viridis() +
  coord_fixed()
out <- levelplot(subset(b, 1:2))

Animate

library(gganimate); library(gifski)
anim <- gplot(subset(b, 1:10)) +
  geom_tile(aes(fill = value)) +
  scale_fill_viridis() +
  coord_fixed() +
  transition_manual(variable) +
  ggtitle("Fuel Moisture Content at {current_frame}") +
  xlab("Longitude") + ylab("Latitude")
animate(anim, duration = 5)

Export an animation

exportdir <- tempdir()
anim_save(filename = "testanimsave.gif", path = exportdir, renderer = gifski_renderer)

Spatial Context Information

I think Open Street Maps can not be downloaded direclty using an R package (see https://github.com/dkahle/ggmap/issues/117). Google maps requires an API key to a google account and possibly for said google account to have a credit card (see help for register_google in package ggmap. Stamen Maps are a possibiity to download using ggmap, however it appears ggmap creates raster objects for these.

Data Supplied by Geoscience Australia

A list of the relevant data is here. Selecting GEODATA TOPO 250K Series 3 opens a new page, defaulting to esri Personal database. One of the 'related' links on the side is for the same data in shape file format. That webpage currently has address [https://ecat.ga.gov.au/geonetwork/srv/eng/catalog.search#/metadata/64058].

According to metadata documents with this data, the projection is latitude and longitude, and the datum is GDA94. The datum GDA94 uses GRS80 as the ellipsoid (Geocentric Datum of Australia 2020 Technical Manual available from here. I'm not sure what other things the datum specifies beyond the ellipsoid. The CRS information extracted from the shape files agrees - they give a proj4string of "+proj=longlat +ellps=GRS80 +no_defs".

railways <- st_read("C:/UserData/hingeek/GA_TOPO_250K_Series_3_shp/Vector_data/Transport/railways.shp")
railways <- railways %>%
  st_transform(st_crs(b)) %>%
  st_crop(b) 

#loading roads below takes a minute or two
roads <- st_read("C:/UserData/hingeek/GA_TOPO_250K_Series_3_shp/Vector_data/Transport/roads.shp")
roads <- roads %>%
  st_transform(st_crs(b)) %>%
  st_crop(b) 
roadsplt <- roads %>%
  ggplot() + 
  geom_sf(aes(col = CLASS))

waterareas <- st_read("C:/UserData/hingeek/GA_TOPO_250K_Series_3_shp/Vector_data/Hydrography/watercourseareas.shp")
waterareas <- waterareas %>%
  st_transform(st_crs(b)) %>%
  st_crop(b) 
waplt <- waterareas %>%
  ggplot() + 
  geom_sf(aes(col = HIERARCHY))

waterlines <- st_read("C:/UserData/hingeek/GA_TOPO_250K_Series_3_shp/Vector_data/Hydrography/watercourselines.shp")
waterlines <- waterlines %>%
  st_transform(st_crs(b)) %>%
  st_crop(b) 
wlplt <- waterlines %>%
  dplyr::filter(HIERARCHY == "Major") %>%
  ggplot() + 
  geom_sf()

builtupareas <- st_read("C:/UserData/hingeek/GA_TOPO_250K_Series_3_shp/Vector_data/Habitation/builtupareas.shp")
builtupareas <- cbind(builtupareas, cen = builtupareas %>% st_centroid() %>% st_coordinates)
builtupareas <- builtupareas %>%
  st_transform(st_crs(b)) %>%
  st_crop(b)

ggplot() +
  geom_sf(data = waterlines %>% dplyr::filter(HIERARCHY == "Major"), stat = "sf", col = "green") + 
  geom_sf(data = roads %>% dplyr::filter(CLASS %in% c("Principal Road", "Dual Carriageway")), inherit.aes = FALSE, stat = "sf", col = "red", lwd = 1) + 
  geom_sf(data = builtupareas %>% dplyr::filter(SHAPE_Area > quantile(builtupareas$SHAPE_Area, 0.8)),
          inherit.aes = FALSE, stat = "sf", fill = "black", col = "black", lwd = 1) +
  coord_sf() +
  geom_text_repel(aes(x = cen.X, y = cen.Y, label = NAME),
        data = builtupareas %>% dplyr::filter(SHAPE_Area > quantile(builtupareas$SHAPE_Area, 0.8)),
        inherit.aes = FALSE, 
        nudge_y = 0.1,
        col = "black",
        size = 3) +
  labs(x = "", y ="")

gplot(subset(b, 1)) +
  geom_tile(aes(fill = value)) +
  scale_fill_viridis() +
  #do not inherit aes in following because don't want to a fill colour given by value
  geom_sf(data = waterareas, inherit.aes = FALSE, stat = "sf", col = "green") + 
  geom_sf(data = waterlines %>% dplyr::filter(HIERARCHY == "Major"), inherit.aes = FALSE, stat = "sf", col = "green") + 
  geom_sf(data = roads %>% dplyr::filter(CLASS == "Principal Road"), inherit.aes = FALSE, stat = "sf", col = "black", lwd = 1) + 
  coord_sf()

Combine into an animation

gplot(subset(b, 1:5)) +
  geom_tile(aes(fill = value)) +
  scale_fill_viridis() +
  #do not inherit aes in following because don't want to a fill colour given by value
  geom_sf(data = waterareas, inherit.aes = FALSE, stat = "sf", col = "green") + 
  geom_sf(data = waterlines %>% dplyr::filter(HIERARCHY == "Major"), inherit.aes = FALSE, stat = "sf", col = "green") + 
  geom_sf(data = roads %>% dplyr::filter(CLASS == "Principal Road"), inherit.aes = FALSE, stat = "sf", col = "black", lwd = 1) + 
  coord_sf() +
  transition_manual(variable) +
  ggtitle("Fuel Moisture Content at {current_frame}") +
  xlab("Longitude") + ylab("Latitude")

SWS Points and GA Context

swspoints <- sws_sites_2_sf(readRDS("./private/data/clean/sws_sites.rds"))

gplot(subset(b, 1:5)) +
  geom_tile(aes(fill = value)) +
  scale_fill_viridis() +
  #do not inherit aes in following because don't want to a fill colour given by value
  geom_sf(data = waterareas, inherit.aes = FALSE, stat = "sf", col = "green") + 
  geom_sf(data = waterlines %>% dplyr::filter(HIERARCHY == "Major"), inherit.aes = FALSE, stat = "sf", col = "green") + 
  geom_sf(data = roads %>% dplyr::filter(CLASS == "Principal Road"), inherit.aes = FALSE, stat = "sf", col = "black", lwd = 1) + 
  geom_sf(data = swspoints, inherit.aes = FALSE, stat = "sf", col = "red") + 
  coord_sf() +
  transition_manual(variable) +
  ggtitle("Fuel Moisture Content at {current_frame}") +
  xlab("Longitude") + ylab("Latitude")

Use Pre-Subsetted GA Data

majorfeatures <- readRDS("./private/data/GA_principalroads_majorrivers_railsways.rds")  %>% 
  st_transform(st_crs(b)) %>%
  st_crop(b)
builtupareas <- readRDS("./private/data/GA_builtupareas.rds")  %>% 
  st_transform(st_crs(b)) %>%
  st_crop(b)

gplot(subset(b, 1:5)) +
  geom_tile(aes(fill = value)) +
  scale_fill_viridis() +
  #do not inherit aes in following because don't want to a fill colour given by value
  geom_sf(aes(col = FEATTYPE), data = majorfeatures, inherit.aes = FALSE, stat = "sf") +
  geom_sf(data = swspoints, inherit.aes = FALSE, stat = "sf", col = "red") + 
  geom_sf(data = builtupareas %>% dplyr::filter(SHAPE_Area > quantile(builtupareas$SHAPE_Area, 0.8)),
          inherit.aes = FALSE, stat = "sf", fill = "black", col = "black", lwd = 1) +
  coord_sf() +
  geom_text_repel(aes(x = cen.X, y = cen.Y, label = NAME),
        data = builtupareas %>% dplyr::filter(SHAPE_Area > quantile(builtupareas$SHAPE_Area, 0.8)),
        inherit.aes = FALSE, 
        nudge_y = 0.1,
        col = "black",
        size = 3) +
  transition_manual(variable) +
  ggtitle("Fuel Moisture Content at {current_frame}") +
  xlab("Longitude") + ylab("Latitude")

And to avoid use of aes for the context (if aes needed for something else):

contextlyrs <- list(
  geom_sf(data = majorfeatures %>% dplyr::filter(FEATTYPE %in% c("Major Road")),
          inherit.aes = FALSE, col = "grey"), 
  # geom_sf(data = majorfeatures %>% dplyr::filter(FEATTYPE %in% c("Major Watercourse")),
  #         inherit.aes = FALSE, lty = "dotted", col = "grey"),
  geom_sf(data = builtupareas %>% dplyr::filter(SHAPE_Area > quantile(builtupareas$SHAPE_Area, 0.8)),
          inherit.aes = FALSE, fill = "grey", col = "grey", lwd = 1),
  geom_text_repel(aes(x = cen.X, y = cen.Y,
                      label = NAME), #to title case: tools::toTitleCase(tolower(as.character(NAME)))),
                  data = builtupareas %>% dplyr::filter(SHAPE_Area > quantile(builtupareas$SHAPE_Area, 0.8)),
                  inherit.aes = FALSE,
                  nudge_y = 0.1,
                  col = "grey",
                  size = 3)
)

gplot(subset(b, 1)) +
  geom_tile(aes(fill = value)) +
  scale_fill_viridis() +
  contextlyrs +
  geom_sf(data = swspoints, inherit.aes = FALSE, stat = "sf", col = "red") + 
  coord_sf() +
  ggtitle("Fuel Moisture Content") +
  xlab("Longitude") + ylab("Latitude")

With Leaflet and NSW Aerial Imagery

library(leaflet)
your.map <- leaflet(swspoints) %>% 
  addTiles() %>%
  ## NSW aerial photography maps aren't useful for this region
  addProviderTiles("Esri.WorldImagery") %>%
  addWMSTiles("https://maps.six.nsw.gov.au/arcgis/services/public/NSW_Imagery/MapServer/WMSServer",
              options = WMSTileOptions(format = "image/png", transparent = F),
              layers = "0") %>%
  addCircleMarkers() %>%
  addScaleBar() %>%
  addMeasure()