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.
Vicmap is the Victorian Government’s catalogue of spatial datasets. The catalogue features 474 datasets across land, property, infrastructure and environment and is the most authoritative suite of spatial data in Victoria.
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.
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.
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.
Note: A previous version of this tutorial used swooping birds as the example, but due to the migration of platforms that data was not available at time of writing
To demonstrate VicmapR’s features, we will download and explore a dataset from the Vicmap catalogue. The data is a spatial dataset of Melbourne’s tram network.
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.
| Name | Title | Abstract | metadataID |
|---|---|---|---|
| open-data-platform:ptv_metro_tram_route | ptv_metro_tram_route | This layer depicts a spatial object (polyline) representing metropolitan tram routes. Each polylines represents a unique route variation. Each route has attributes that describe route, trip headsign (direction), route length, first/last stop, number of stops and operator name. The data has been generated from the PTV GTFS data with extra attributes from PTV’s TransNet database. This dataset supersedes "PTV_TRAM_ROUTE". This dataset was first loaded into the VSDL in March 2018 and will be updated approximately quarterly. | be15bc57-a614-5df1-a2a7-2e8d39fe6d1b |
| open-data-platform:ptv_metro_tram_stop | ptv_metro_tram_stop | This layer depicts spatial objects (points) representing metropolitan tram stops. Each tram stop has attributes that describe StopID, StopName, Lat/Long, Ticket Zone and Routes Using Stop. The data has been generated from the PTV GTFS data with extra attributes from PTV’s TransNet database. This dataset supersedes "PTV_TRAM_STOP" (which included metro and regional stops). This dataset was first loaded into the VSDL in March 2018 and will be updated approximately quarterly. | 0058eda4-7b1e-5681-9758-7adbda8481b1 |
| open-data-platform:ptv_tram_track_centreline | ptv_tram_track_centreline | Tram Track Centreline depicts a spatial object (polyline) representing the track centreline (the centre of the two rails). | c83faa88-cc6b-57ab-a86c-969c324d8adf |
| open-data-platform:recweb_historic_relic | recweb_historic_relic |
Recreation historic relic dataset describes historic relics suitable for public visitation (such as timber tramways, old sawmill sites etc) that are promoted for visitation. The recreation historic relics within State Forest have been captured and recorded with a Trimble Pro XR GPS and are actively promoted to the public and maintained by the Department of Environment, Land, Water and Planning.
|
f820e396-84ac-598c-9c73-f0305a028783 |
| open-data-platform:tr_rail | tr_rail | This layer is part of Vicmap Transport and contains line and point features delineating railway features. Includes; Railway Exits, Railway Yards, Railway Bridges, Railway Tunnels, Railways & Tramways. | 0419297b-9552-5797-89f0-e7f8ca251e7d |
| open-data-platform:vic_tourist_rail_corridor | vic_tourist_rail_corridor |
This polyline dataset represents the centreline of active
tourist/heritage railway corridors in Victoria.
The dataset inclues the Bellarine Penninsula Tourist Railway, Daylesford Spa Country Tourist Railway, Mornington Tourist Railway, Puffing Billy Tourist Railway, South Gippsland Tourist Railway, Victorian Goldfields Railway, Walhalla Goldfields Tourist Railway and Yarra Valley Tourist Railway. The railways are run by non-profit groups who are mainly volunteers and rely on strong community involvement and support. The dataset does not include several miniture railways and tramways also found in Victoria. |
6827bc20-8a84-5175-b0f8-4bf5151a022e |
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")
#> • Using collect() on this object will return 115 features
#> • and 14 fields
#> • At most six rows of the record are printed here
#> ────────────────────────────────────────────────────────────────────────────────
#> Simple feature collection with 6 features and 13 fields
#> Geometry type: LINESTRING
#> Dimension: XY
#> Bounding box: xmin: 144.942 ymin: -37.87267 xmax: 145.0286 ymax: -37.74284
#> Geodetic CRS: GDA94
#> # A tibble: 6 × 14
#> id shape_id route_id route_short_name trip_headsign route_km
#> <chr> <chr> <chr> <chr> <chr> <dbl>
#> 1 ptv_metro_tram_rout… 3-12-A-… 3-12-A-… 12 Victoria Gar… 5.59
#> 2 ptv_metro_tram_rout… 3-12-A-… 3-12-A-… 12 Victoria Gar… 11.4
#> 3 ptv_metro_tram_rout… 3-12-A-… 3-12-A-… 12 Victoria Gar… 6.62
#> 4 ptv_metro_tram_rout… 3-16-A-… 3-16-A-… 16 Kew to Melbo… 14.4
#> 5 ptv_metro_tram_rout… 3-11-A-… 3-11-A-… 11 West Preston… 10.2
#> 6 ptv_metro_tram_rout… 3-1-A-m… 3-1-A-m… 1 South Melbou… 1.21
#> # ℹ 8 more variables: route_long_name <chr>, first_stop_name <chr>,
#> # last_stop_name <chr>, first_stop_id <chr>, last_stop_id <chr>,
#> # num_of_stops <int>, operator_name <chr>, geometry <LINESTRING [°]>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()| Metadata Name | Descriptions |
|---|---|
| Resource Name | PTV_METRO_TRAM_ROUTE |
| Title | PTV Metro Tram Routes |
| Anzlic ID | ANZVI0803005849 |
| Custodian | NA |
| Owner | NA |
| Jurisdiction | Victoria |
| Abstract | This layer depicts a spatial object (polyline) representing metropolitan tram routes. Each polylines represents a unique route variation. Each route has attributes that describe route, trip headsign (direction), route length, first/last stop, number of stops and operator name. The data has been generated from the PTV GTFS data with extra attributes from PTV’s TransNet database. This dataset supersedes “PTV_TRAM_ROUTE”. This dataset was first loaded into the VSDL in March 2018 and will be updated approximately quarterly. |
| Search Words | Transportation |
| Purpose | This data is intended to be used by any application which needs to show the location of metropolitan tram routes in Victoria. |
| Geographic Extent Polygon | NA |
| Geographic Bounding Box | NA |
| Beginning to Ending Date | 2018-03-01 - now |
| Maintainence and Update Frequency | Quarterly |
| Stored Data Format | NA |
| Available Format(s) Types | DIGITAL |
| Positional Accuracy | This data was generated from GTFS data. Some spatial inaccurancies may exist. |
| Attribute Accuracy | Not Applicable |
| Logical Consistency | NA |
| Data Source |
Dataset Source: This data was generated using a MapBasic process from PTV GTFS data and an export from TranNet. Dataset Originality: Derived |
| Contact Organisation | NA |
| Contact Position | NA |
| Address | NA |
| Telephone | NA |
| Facsimile | NA |
| Email Address | NA |
| Metadata Date | 2022-02-16 |
| Additional Metadata |
Related Documents: None Not specified |
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
#> Simple feature collection with 10 features and 13 fields
#> Geometry type: LINESTRING
#> Dimension: XY
#> Bounding box: xmin: 144.942 ymin: -37.87267 xmax: 145.0286 ymax: -37.72939
#> Geodetic CRS: GDA94
#> # A tibble: 10 × 14
#> id shape_id route_id route_short_name trip_headsign route_km
#> <chr> <chr> <chr> <chr> <chr> <dbl>
#> 1 ptv_metro_tram_rou… 3-12-A-… 3-12-A-… 12 Victoria Gar… 5.59
#> 2 ptv_metro_tram_rou… 3-12-A-… 3-12-A-… 12 Victoria Gar… 11.4
#> 3 ptv_metro_tram_rou… 3-12-A-… 3-12-A-… 12 Victoria Gar… 6.62
#> 4 ptv_metro_tram_rou… 3-16-A-… 3-16-A-… 16 Kew to Melbo… 14.4
#> 5 ptv_metro_tram_rou… 3-11-A-… 3-11-A-… 11 West Preston… 10.2
#> 6 ptv_metro_tram_rou… 3-1-A-m… 3-1-A-m… 1 South Melbou… 1.21
#> 7 ptv_metro_tram_rou… 3-11-A-… 3-11-A-… 11 Victoria Har… 13.3
#> 8 ptv_metro_tram_rou… 3-12-A-… 3-12-A-… 12 St Kilda (Fi… 6.61
#> 9 ptv_metro_tram_rou… 3-12-A-… 3-12-A-… 12 St Kilda (Fi… 11.4
#> 10 ptv_metro_tram_rou… 3-12-A-… 3-12-A-… 12 St Kilda (Fi… 5.62
#> # ℹ 8 more variables: route_long_name <chr>, first_stop_name <chr>,
#> # last_stop_name <chr>, first_stop_id <chr>, last_stop_id <chr>,
#> # num_of_stops <int>, operator_name <chr>, geometry <LINESTRING [°]>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()| route_id | route_short_name | trip_headsign | route_km | operator_name | geometry |
|---|---|---|---|---|---|
| 3-12-A-mjp-1 | 12 | Victoria Gardens to St Kilda (Fitzroy St) | 11.39 | Yarra Trams | LINESTRING (145.0101 -37.81… |
| 3-12-A-mjp-1 | 12 | Victoria Gardens to St Kilda (Fitzroy St) | 6.62 | Yarra Trams | LINESTRING (145.0101 -37.81… |
| 3-16-A-mjp-1 | 16 | Kew to Melbourne University via St Kilda Beach | 14.38 | Yarra Trams | LINESTRING (145.0286 -37.86… |
| 3-11-A-mjp-1 | 11 | West Preston to Victoria Harbour Docklands | 10.20 | Yarra Trams | LINESTRING (144.995 -37.752… |
| 3-11-A-mjp-1 | 11 | Victoria Harbour Docklands to West Preston | 13.31 | Yarra Trams | LINESTRING (144.942 -37.820… |
By filtering the Vicmap_promise before
collecting, we have downloaded only a subset of rows of
data.
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.
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)
)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!
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 tram route 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 = "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)
#> [1] TRUE
# Intersection
darebin_trams <- sf::st_intersection(tram_route, darebin)
mapview(darebin_trams) + mapview(darebin, alpha.regions = 0.3, col.regions= "green")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 tram route data. Don’t forget, we need to convert the crs to 4283 to work with our tram route data
st_crs(home_radius) #crs is UTM Zone 55
#> Coordinate Reference System:
#> User input: +proj=utm +zone=55
#> wkt:
#> PROJCRS["unknown",
#> BASEGEOGCRS["unknown",
#> DATUM["World Geodetic System 1984",
#> ELLIPSOID["WGS 84",6378137,298.257223563,
#> LENGTHUNIT["metre",1]],
#> ID["EPSG",6326]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["degree",0.0174532925199433],
#> ID["EPSG",8901]]],
#> CONVERSION["UTM zone 55N",
#> METHOD["Transverse Mercator",
#> ID["EPSG",9807]],
#> PARAMETER["Latitude of natural origin",0,
#> ANGLEUNIT["degree",0.0174532925199433],
#> ID["EPSG",8801]],
#> PARAMETER["Longitude of natural origin",147,
#> ANGLEUNIT["degree",0.0174532925199433],
#> ID["EPSG",8802]],
#> PARAMETER["Scale factor at natural origin",0.9996,
#> SCALEUNIT["unity",1],
#> ID["EPSG",8805]],
#> PARAMETER["False easting",500000,
#> LENGTHUNIT["metre",1],
#> ID["EPSG",8806]],
#> PARAMETER["False northing",0,
#> LENGTHUNIT["metre",1],
#> ID["EPSG",8807]],
#> ID["EPSG",16055]],
#> CS[Cartesian,2],
#> AXIS["(E)",east,
#> ORDER[1],
#> LENGTHUNIT["metre",1,
#> ID["EPSG",9001]]],
#> AXIS["(N)",north,
#> ORDER[2],
#> LENGTHUNIT["metre",1,
#> ID["EPSG",9001]]]]
home_radius <- st_transform(home_radius, crs = 4283)
st_crs(home_radius)
#> Coordinate Reference System:
#> User input: EPSG:4283
#> wkt:
#> GEOGCRS["GDA94",
#> DATUM["Geocentric Datum of Australia 1994",
#> ELLIPSOID["GRS 1980",6378137,298.257222101,
#> LENGTHUNIT["metre",1]]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["degree",0.0174532925199433]],
#> CS[ellipsoidal,2],
#> AXIS["geodetic latitude (Lat)",north,
#> ORDER[1],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> AXIS["geodetic longitude (Lon)",east,
#> ORDER[2],
#> ANGLEUNIT["degree",0.0174532925199433]],
#> USAGE[
#> SCOPE["Horizontal component of 3D system."],
#> AREA["Australia including Lord Howe Island, Macquarie Island, Ashmore and Cartier Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All onshore and offshore."],
#> BBOX[-60.55,93.41,-8.47,173.34]],
#> ID["EPSG",4283]]
tram_route_10k <- st_intersection(tram_route, home_radius)
mapview(tram_route_10k) + mapview(home_radius, alpha.regions = 0.3, col.regions = "green")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 get caught in tram tracks.
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)
#> Reading layer `cycle_route' from data source
#> `/Users/runner/work/_temp/Library/VicmapR/shapes/cycle_route.geojson'
#> using driver `GeoJSON'
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: LINESTRING
#> Dimension: XY
#> Bounding box: xmin: 144.9627 ymin: -37.8303 xmax: 145.1933 ymax: -37.7021
#> Geodetic CRS: WGS 84
#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)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 VicmapR 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 tram routes only in
Darebin we can filter with INTERSECTS().
# trams in Darebin
darebin_trams <- vicmap_query(layer = "open-data-platform:ptv_metro_tram_route") %>%
filter(INTERSECTS(darebin)) %>%
collect()
mapview(darebin_trams) + mapview(darebin, alpha.regions = 0.3, col.regions = "orange")And finally, the tram routes on our cycling route. To find tram
tracks within 30m 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)
#> Reading layer `cycle_route' from data source
#> `/Users/runner/work/_temp/Library/VicmapR/shapes/cycle_route.geojson'
#> using driver `GeoJSON'
#> Simple feature collection with 1 feature and 0 fields
#> Geometry type: LINESTRING
#> Dimension: XY
#> Bounding box: xmin: 144.9627 ymin: -37.8303 xmax: 145.1933 ymax: -37.7021
#> Geodetic CRS: WGS 84
# Condense line object
route_poly <- route_line %>%
sf::st_transform(3111) %>%
sf::st_buffer(30) %>%
sf::st_cast("POLYGON") %>%
sf::st_transform(4283)
route_intersection <- vicmap_query(layer = "open-data-platform:ptv_metro_tram_route") %>%
filter(INTERSECTS(route_poly)) %>%
collect() %>%
st_intersection(route_poly)
mapview::mapview(route_poly) + mapview::mapview(route_intersection, color = "red", lwd = 5)