Vicmap for Beginners

2022-06-10

Who is this tutorial for?

The purpose of this tutorial is:

The end goal of this tutorial is for R users with minimal spatial experience to be able to search for, download and visualise Vicmap data. While this tutorial is aimed at spatial beginners, it can still be used as a reference for more experienced users.

What is Vicmap?

Vicmap is the Victorian Government’s catalogue of spatial datasets. The catalogue features 651 datasets across land, property, infrastructure and environment and is the most authoritative suite of spatial data in Victoria.

Accessing Vicmap datasets

This data catalogue is freely available to the public and can be accessed via several methods. For users of R, the fastest way to access up to date Vicmap datasets is to utilise the Web Feature Service (WFS). WFS is a standardised interface to request geographic information, regardless of the platform on which it is stored.

WFS requires a URL that contains the instructions for the query and is written in WFS specific terminology. This of course requires understanding of the WFS terminology and the time-consuming process of manually building the URL string.

Enter, VicmapR.

WFS, without the fuss!

VicmapR simplifies WFS queries by taking user supplied keywords and automatically building the correctly formatted URL string.
The functions in VicmapR therefore provide a much faster and more robust method for acquiring Vicmap data.

Lazy evaluation

VicmapR uses lazy evaluation, in which we query but do not collect the data from the Vicmap database. This is called making a ‘Vicmap Promise’. A Vicmap Promise contains the dataset information but not the data itself. This promise can be filtered for a subset of the data before collecting. Using lazy evaluation we can collect the desired subset, instead of the entire dataset.

Worked Example - Swooping Birds

To demonstrate VicmapR’s features, we will download and explore a dataset from the Vicmap catalogue.

Australian spring is famous for the presence of swooping birds. Certain species of birds swoop during the breeding season to protect their nest from potential predators. Swooping birds in Victoria include the masked lapwing (plover), butcherbird, magpie-lark, noisy miner and perhaps most famously, the magpie. The Vicmap catalogue has a dataset of ‘swoop’ locations, recorded by members of the public.

Searching for data

First, we need to search for the dataset. We can view all layers, or do a keyword search. Use ignore.case = TRUE for a case-insensitive search.

#Don't forget to load VicmapR library
library(VicmapR) 

#All layers
all_layers <- listLayers() 

#Case insensitive keyword search
search <- listLayers(pattern = "swoop", ignore.case = T)

Viewing our search results shows us the name of the dataset and its description.

Once we have the name of the Vicmap dataset we would like to download, in this case ‘datavic:FLORAFUANA1_SWOOPING_BIRD’, we can query the Vicmap database.

If we are still not sure from the name and description that this is the right dataset, we can look at the Vicmap Promise to get a snippet of the dataset contents.

vicmap_query(layer = "datavic:FLORAFUANA1_SWOOPING_BIRD")

We can see that the dataset contains, swoop date, location and species. At a glance this data should be adequate to start analysing the spatial distribution of bird swoops.

Downloading the data

The summary also shows that there are rows in this dataset. This isn’t a very large dataset, but if we are using this data in an interactive report or application, waiting to download the full dataset may be impractical.

Using lazy evaluation, we can use pipes to filter and subset the data, so that we only collect the desired subset. For example, let’s just look at the first 10 rows.

query <- vicmap_query(layer = "datavic:FLORAFUANA1_SWOOPING_BIRD") %>% 
  head(10) %>% collect()

query

Subsetting data

Now let’s only look at swoops that occurred in 2021, where STATUS is HISTORY2021. For our analysis, we will look only at the swoop date, species, number of birds, comments and geometry. Note that Vicmap datasets will retain id and OBJECTID columns even if not selected in the Vicmap promise. As this a simple features (sf) dataset, the geometry corresponding to the features (in this case, swoops) does not need to be selected as it is always attached. The geometry column can only be removed intentionally, for example, by using the st_drop_geometry() function in the sf package.

library(dplyr)

swoops_2021 <- vicmap_query(layer = "datavic:FLORAFUANA1_SWOOPING_BIRD") %>% 
  filter(SWOOP_DATE > "2020-12-31z") %>% 
  filter(STATUS == "HISTORY2021") %>% 
  select(SWOOP_DATE, NO_OF_BIRDS, SPECIES, COMMENTS) %>%
  collect()

swoops_2021 %>% head(5) %>% 
  select(-id) %>% #condense table for viewing
  kable() %>% 
  kable_styling()

By filtering the Vicmap promise before collecting, we have downloaded only a subset of rows of data.

Quick visualisation of spatial data with the mapview package

Often the first thing we want to do is see the data on a map. The simplest and quickest way to visualise spatial data is to use the mapview package:

library(mapview)
mapview(swoops_2021)

While the output is basic, it is very simple to achieve and points can be clicked on to obtain the details from the other fields.

Visualisation with leaflet

For more customisation, you might want to consider plotting your map in leaflet. Leaflet is JavaScript library that allows users to produce interactive maps that work efficiently across most desktop and mobile platforms. Using the leaflet package in R, we have the flexibility of plotting more complex map products.

We will start by plotting the most basic map in leaflet. Instead of using a single function like mapview(), in leaflet we need to add each component to the map. At a minimum we specify the leaflet object, add a basemap and add the markers. Note: differing to mapview, leaflet does not include marker labels by default, these must be added with the popup parameter.

library(leaflet)

swoops_2021 %>%
  leaflet() %>% #plot a leaflet object
  addTiles() %>% #add a basemap (default is OpenStreetMap)
  addMarkers(popup = swoops_2021$COMMENTS) #add geom points the marker label (automatically obtained from geom column)

There are a lot of points in this dataset and leaflet’s default marker appears very large and crowded. addCircleMarkers() is a quick and easy way to customise markers to suit your data.

swoops_2021 %>%
  leaflet() %>%
  addTiles() %>%
  addCircleMarkers(
    radius = 6,
    color = "green",
    stroke = FALSE, #removes outline
    fillOpacity = 0.6,
    popup = swoops_2021$COMMENTS
  )

We can also easily conditionally format our markers and basemaps in leaflet. Let’s:


# Create a palette
pal <- colorFactor("Accent", levels = swoops_2021$SPECIES) #define a colour palette

swoops_2021 %>%
  leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>% #add third party base map
  addCircleMarkers(
    radius = ~NO_OF_BIRDS,
    color = ~pal(SPECIES),
    stroke = FALSE,
    fillOpacity = 0.7,
    popup = paste0("<b>COMMENT: </b>", swoops_2021$COMMENTS, "<br>", #format the popup with html tags
                   "<b>NUM. BIRDS: </b>", swoops_2021$NO_OF_BIRDS, "<br>",
                   "<b>SPECIES: </b>", swoops_2021$SPECIES)
  )

Spatial filtering

Let’s say we want to know about the swoop spots only in our LGA (e.g. Darebin). We can apply a spatial filter to our dataset so that only the swoops in our suburb are shown. Where do we get the spatial data for our LGA? From the Vicmap catalogue of course!

The sf package

The sf package for R provides tools for dealing with ‘simple feature’ datasets, i.e. a data frame or tibble with a geometry list-column. The sf package allows us to manupulate spatial data and perform spatial operations. In this case, we want to find which points in our swoops data intersect with the polygon for our suburb polygon.

To do this we can use sf::st_intersection(). Before performing any spatial operations on a dataset, you want to make sure the coordinate reference system (crs) is correct and if there are multiple datasets, that the crs is the same for each dataset. You can check the crs system with sf::st_crs().

# Suburb polygon for Darebin
darebin <- vicmap_query(layer = "datavic:VMADMIN_LGA_POLYGON") %>%
  select(LGA_NAME) %>% 
  filter(LGA_NAME == "DAREBIN") %>%
  collect()


library(sf)

# Check crs for each spatial dataset
sf::st_crs(darebin) == sf::st_crs(swoops_2021)

# Intersection
darebin_swoops <- sf::st_intersection(swoops_2021, darebin)
mapview(darebin_swoops) + mapview(darebin, alpha.regions = 0.3, col.regions= "green")

Filtering within a radius

We probably regularly travel a bit further afield though, so lets look at a radius from our home address. We can create the radius geometry with st_buffer() but first need to define the coordinates for the centre of the radius (home). Don’t forget to specify the crs as 4326 as we are dealing with a latitude and longitude coordinate system.

lat <- -37.75215888969604
lon <- 145.02927170548745

home <- st_sfc(st_point(c(lon, lat)), crs = 4326)

We have a problem though, latitude and longitude units are degrees and st_buffer assumes units of meters. To convert from degrees to metres we need to project our latitude and longitude, which are represented as a point on a curved surface, onto a flat plane. To do this, we use a projection standard chosen specifically for the zone containing the coordinates. The standard is called the EPSG code and for our coordinates is zone 55.

# Project coordinates
home_utm <- st_transform(home, "+proj=utm +zone=55")

# Get buffer
home_radius <- sf::st_buffer(home_utm, dist = 10000)

mapview(home_utm) + mapview(home_radius)

Now we can use this radius to filter our swoops data. Don’t forget, we need to convert the crs to 4283 to work with our swoops data

st_crs(home_radius) #crs is UTM Zone 55

home_radius <- st_transform(home_radius, crs = 4283)
st_crs(home_radius)

swoops_10k <- st_intersection(swoops_2021, home_radius)
mapview(swoops_10k) + mapview(home_radius, alpha.regions = 0.3, col.regions = "green")

Finding nearby points

Let’s say we feel like a bike ride. We’ve picked a route but want to know which areas to avoid so as not to disturb too many nesting birds. This is a more difficult problem as the ride travels through several LGAs and we only really want to know about reported swoops that are very close to the route.

We can solve this problem in a similar way, using st_intersection().


# Load bike route
route_line <- sf::st_read(system.file("shapes/cycle_route.geojson", package="VicmapR"), quiet = F) %>% 
  sf::st_transform(4283)

#Add buffer to route - we need to convert to m to do this (same as before)
route_m <- st_transform(route_line, "+proj=utm +zone=55") %>% st_buffer(300) #300m buffer
swoops_m <- st_transform(swoops_2021, "+proj=utm +zone=55") # convert to same crs for st_intersection

swoops_en_route <- st_intersection(route_m,swoops_m)

# Recorded swoops within 300m of the route
mapview(swoops_en_route) + mapview(route_line)

Using the VicmapR Geometric Filters

The VicmapR package offers tools for geometric filtering that avoid the need for sf package. The WFS Geoserver on which Vicmap is based supports several geometric filters - see the full list here.

These geometric filters allow us to perform basic spatial manipulation while retaining the benefits of lazy evaluation. So as before, if we want to filter our data before downloading, we can do so with more complex spatial operations. This is particularly useful when presenting Vicmap data in an interactive map, as it reduces the time for downloading data between re-rendering the map.

Let’s revisit our examples above. To get the swoops only in Darebin we can filter with INTERSECTS().


# Suburb polygon for Darebin
darebin <- vicmap_query(layer = "datavic:VMADMIN_LGA_POLYGON") %>%
  select(LGA_NAME) %>% 
  filter(LGA_NAME == "DAREBIN") %>%
  collect()

# Swoops in Darebin
darebin_2021 <- vicmap_query(layer = "datavic:FLORAFUANA1_SWOOPING_BIRD") %>%
  filter(SWOOP_DATE > "2020-12-31z") %>%
  filter(STATUS == "HISTORY2021") %>%
  select(SWOOP_DATE, NO_OF_BIRDS, SPECIES, COMMENTS) %>%
  filter(INTERSECTS(darebin)) %>%
  collect()

mapview(darebin_2021) + mapview(darebin, alpha.regions = 0.3, col.regions = "orange")

Now let’s try the 10km radius.


# Swoops in 10k
radius_2021 <- vicmap_query(layer = "datavic:FLORAFUANA1_SWOOPING_BIRD") %>%
  filter(SWOOP_DATE > "2020-12-31z") %>%
  filter(STATUS == "HISTORY2021") %>%
  select(SWOOP_DATE, NO_OF_BIRDS, SPECIES, COMMENTS) %>%
  filter(INTERSECTS(home_radius)) %>%
  collect()

mapview(radius_2021) + mapview(home_radius, alpha.regions = 0.3, col.regions = "orange")

And finally, the swoops on our cycling route. To find swoops within 300m we use the INTERSECTS() function and the cycling route with an added 300m buffer. Note: VicmapR geometric filters will be simplified, thus if precise intersections are desired, it is advised to do a cleanup of the filtered data once collected.


route_line <- sf::st_read(system.file("shapes/cycle_route.geojson", package="VicmapR"), quiet = F) %>% 
  sf::st_transform(4283)

# Condense line object
route_poly <- route_line %>% 
  sf::st_transform(3111) %>% 
  sf::st_buffer(300) %>% 
  sf::st_cast("POLYGON") %>% 
  sf::st_transform(4283)

swoop <- VicmapR::vicmap_query("datavic:FLORAFUANA1_SWOOPING_BIRD") %>%
  filter(SWOOP_DATE > "2020-12-31z") %>%
  filter(STATUS == "HISTORY2021") %>%
  select(SWOOP_DATE, NO_OF_BIRDS, SPECIES, COMMENTS) %>% 
  filter(INTERSECTS(route_poly)) %>% collect()


mapview::mapview(route_poly) + mapview::mapview(swoop)

Debrief

After completing this tutorial you should have a basic understanding of plotting spatial data and performing some basic spatial operations. This tutorial is a quick guide and designed only to get you started. For further learning, see the following resources: