Vicmap for Beginners

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.

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.

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. listLayers() by default now provides a column of ‘Abstract’ and ‘metadataID’. You can remove these columns by setting the abstract argument to FALSE (e.g. listLayers(abstract = FALSE))

#All layers
all_layers <- listLayers() 

# Case insensitive keyword search
search <- listLayers(pattern = "tram", 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 ‘open-data-platform:ptv_metro_tram_route’, 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 = "open-data-platform:ptv_metro_tram_route") 

We can see that the dataset contains the route numbers and names. At a glance this data should be adequate. If you want more extensive information of the data you can use get_metadata() to download a list of (i) metadata about the data, (ii) a data dictionary (work-in-progress) and (iii) a link to the metadata url.

metadata <- vicmap_query(layer = "open-data-platform:ptv_metro_tram_route") %>%
  get_metadata() %>%
  .[[1]]

kbl(metadata) %>%
  kable_styling()

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 = "open-data-platform:ptv_metro_tram_route") %>% 
  head(10) %>% collect()

query

Subsetting data

Now let’s only look at the tram route data. Note that VicmapR datasets will retain id and geometry columns even if not selected in the VicmapR promise. As this a simple features (sf) dataset, the geometry corresponding to the features 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.

tram_select <- vicmap_query(layer = "open-data-platform:ptv_metro_tram_route") %>% 
  filter(operator_name == "Yarra Trams") %>% 
  filter(num_of_stops > 20) %>% 
  select(route_id, route_short_name, trip_headsign, operator_name, route_km) %>%
  collect()

tram_select %>% 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)
mapviewOptions(fgb = FALSE)

# Because the dataset is not big we can download all the data as 'tram_route'
tram_route <- vicmap_query(layer = "open-data-platform:ptv_metro_tram_route") %>% 
  collect() 

# plot data
mapview(tram_route)

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

Visualisation with leaflet

For more customization, 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)

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

tram_route %>%
  leaflet() %>%
  addProviderTiles(providers$CartoDB.Positron) %>% #add third party base map
  addPolylines(
    color = ~pal(trip_headsign),
    weight = 2,
    stroke = 0.5, #removes outline
    fillOpacity = 0.8,
    popup = paste0("<b>Route No.: </b>", tram_route$route_short_name, "<br>", #format the popup with html tags
                   "<b>Route Name: </b>", tram_route$trip_headsign, "<br>",
                   "<b>Distance (km): </b>", tram_route$route_km)
  )

Spatial filtering

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

# Suburb polygon for Darebin
darebin <- vicmap_query(layer = "open-data-platform:lga_polygon") %>%
  select(lga_name) %>% 
  filter(lga_name == "DAREBIN") %>%
  collect() %>%
  st_make_valid() # magic fix for some spatial data

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

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

lat <- -37.75215888969604
lon <- 145.02927170548745

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

# 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)

st_crs(home_radius) #crs is UTM Zone 55

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

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

# 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")
tram_m <- st_transform(tram_route, "+proj=utm +zone=55") # convert to same crs for st_intersection

trams_en_route <- st_intersection(route_m,tram_m)

# Recorded tram tracks within 300m of the route
mapview(trams_en_route, col.regions = "Black") + mapview(route_line)