Basemaps, eh?: A guide to creating Canadian basemaps using {rcanvec}, {rosm} and {prettymapr}

by on November 22, 2015

Install the packages

See if you can install the packages – they’re on CRAN and should be available for most verions of R. Updating R is a pain, but it may be what you have to do to make sure you can install the packages. If you already have the packages you’ll want to do this anyway to make sure you have the latest version of all of them.

```install.packages("prettymapr")
install.packages("rosm")
install.packages("rcanvec")
```

Later on you’ll also need some data, which you can find in the tutorial notes.

Base plotting

Part of being easy to use means these packages work well with the {sp} package, or any other package that uses base plotting in R (e.g. {marmap}, {cartography}, {OpenStreetMap}). Because of this, it’s worth taking a look at how to use the plot() function and a few related methods. First, we’ll setup the data we’ll use in our examples.

```# set up data
xs <- c(-64.35134, -64.36888, -64.35984)
ys <- c(45.06345, 45.08933, 45.09176)
mylabels <- c("Graveyard", "Huggins", "TAN")
```

plot()

The `plot(x, y)` command creates a new plot and plots its arguments as points. Usually {rcanvec} or {rosm} will take care of creating the initial plot, but it’s worth knowing that behind the scenes, this is what’s happening. You can use `xlim=c(X, X)`, and `ylim=c(X, X)` to change the extents of the plot.

```plot(xs, ys)
```

Now that we have our basic plot, we’ll use a few other functions to draw on top of the plot we have just created.

points()

It’s a bit redundant to add points on top of the ones we’ve just created, but often you’ll want to add points to a plot without creating a new one. The syntax is exactly the same as for `plot()`. Use `cex` (or character expantion factor) to make your points bigger, `pch` (or ‘plotting character’; see this Quick-R page for options) to change the symbol, or `col` to change the colour (try “black”, “white”, “#554499”, etc.).

```# initial plot
plot(xs, ys)
points(xs, ys, pch=2)
```

arrows()

This does exactly what you’d think: draw an arrow from `x1`, `y1` to `x2`, `y2`. Use the `length` argument to make the arrows bigger or smaller (in inches).

```# initial plot
plot(xs, ys)
# draw an arrow from "Graveyard" to "Huggins"
arrows(xs[1], ys[1], xs[2], ys[2])
```

lines()

If you’d like to connect your points with lines, the `lines()` function is for you! In addition to the `col` argument, we can also pass `lwd` (line width; try from 0.5 to 2), or lty (line type; you can pick any number from 0 to 6).

```# initial plot
plot(xs, ys)
# add lines from the graveyard to huggins to TAN.
lines(xs, ys)
```

text()

Labels are often useful, and adding them to the map is no problem. Use cex to control font size (try from 0.5 to 2), col to control colour, and adj to offset the label. Generally you will want to use `adj=c(-0.2, 0.5)`, which will offset the label to the right of your points.

```# initial plot
plot(xs, ys)
text(xs, ys, labels=mylabels)
```

polygon()

As advertized, creates a polygon based on the verticies provided. Try border to change the border colour, col to change the fill colour, lty to change the border line type, and lwd to change the border line width.

```# initial plot
plot(xs, ys)
# add a polygon between the three points
polygon(xs, ys, col="green", border="red")
```

locator()

Often it’s desirable to place something on the plot that you don’t necessarily know the coordinates to. To find these coordinates, we can capture mouse clicks on a plot using the locator() function. This will return a data.frame of XY coordinates of the points we clicked on. Remember to hit ESC when done!

```# use locator() to find points (hit ESC when done!)
locator()
```

It’s also possible to capture this to a variable, but you’ll eventually want to somehow hardcode these values into your script (because you won’t want to click on the map every single time you run your script).

More Resources

There are numerous graphical parameters you can pass to these functions, which are documented in great detail at Quick-R and in the R man page for par. The same parameters will be used when you plot spatial data, so it’s worth becoming familiar with how to make points, lines, and polygons the way you’d like them to look.

Plotting Spatial Data

The {rgdal} and {sp} packages provide a powerful engine to render all kinds of geographical data. The data used in this section can be found in the tutorial notes you’ve hopefully already downloaded. First, we’ll have to load the packages.

```#load packages
library(sp)
library(rgdal)
```

Next, we’ll load the data. `readOGR()` from the {rgdal} package lets us load shapefiles, among many other data types. The syntax is a bit tricky, but for shapefiles dsn is the folder that contains your shapefile, and layer is the name of the shapefile (without the “.shp”).

```#use rgdal::readOGR() to read shapefile
#dsn=SOURCE FOLDER, #layer=FILENAME (minus .shp)
altashoreline <- readOGR(dsn="data", layer="alta_shoreline")
```
```## OGR data source with driver: ESRI Shapefile
## Source: "data", layer: "alta_shoreline"
## with 1 features
## It has 23 fields```

Once we’ve loaded the data, we can plot it. Use the graphical parameters such as col, border, and lwd we discussed earlier to make the data look decent. You can also use `axes`, `xlim`, and `ylim` to customize the appearance of the plot and extents, respectively. With this plot command we can also use `add=TRUE` to plot another file on top (instead of creating a new plot).

```# plot data
plot(altashoreline, col="lightblue")
```

Now that we’ve plotted our base data, we can draw on top of it using `lines()`, `polygon()`, `arrows()`, and `text()`. As an example, we’ll load a file containing the locations of my master’s thesis sample locations.

```# plot shoreline
plot(altashoreline, col="lightblue")
# plot Alta Lake and cores
points(altacores\$lon, altacores\$lat, pch=15, cex=0.9)
text(altacores\$lon, altacores\$lat, labels=altacores\$name, adj=c(-0.2, 0.5), cex=0.9)
```

Using {rosm} to plot basemaps

The {rosm} package pulls Bing Maps, Open Street Map, and related maps from the internet, caches them locally, and renders them to provide context to overlying data (your sample sites, etc.). Again, for details, take a look at the {rosm} manual. First we’ll load the packages.

```library(prettymapr)
library(rosm)
```

Step 1: Find your bounding box

The {rosm} package plots based on a bounding box, or an area that you would like to be visible. There’s a few ways to go about doing this, but the easiest way is to visit the Open Street Maps Export page, zoom to your area of interest, and copy/paste the values into `makebbox(northlat, eastlon, southlat, westlon)` from the {prettymapr} package. You can also use `searchbbox("my location name")`, also from the {prettymapr} package, which will query google for an appropriate bounding box. You’ll notice that the bounding box returned by these methods is just a 2×2 matrix, the same as that returned by `bbox()` in the {sp} package.

```altalake <- searchbbox("alta lake, BC")
# or
altalake <- makebbox(50.1232, -122.9574, 50.1035, -123.0042)
```

Make sure you’ve got your bounding box right by trying `osm.plot()` or `bmaps.plot()` with the bounding box as your first argument.

```osm.plot(altalake)
```

```bmaps.plot(altalake)
```

Step 2: Choose your map type and zoom level

{rosm} provides access to a number of map types (and even the ability to load your own if you’re savvy!), but the most common ones you’ll use are `type=osm`, `type="hillshade"`, `type=stamenwatercolor`, and `type=stamenbw` for `osm.plot()` and `type="Aerial"` with `bmaps.plot()`. Look at all of them with `osm.types()` and `bmaps.types()`.

```osm.types()
```
```##  [1] "hikebike"               "hillshade"
##  [3] "hotstyle"               "lovinacycle"
##  [5] "lovinahike"             "mapquestosm"
##  [7] "mapquestsat"            "opencycle"
##  [9] "openpiste"              "osm"
## [11] "osmgrayscale"           "osmtransport"
## [13] "stamenbw"               "stamenwatercolor"
## [15] "thunderforestlandscape" "thunderforestoutdoors"```
```osm.plot(altalake, type="stamenbw")
```

```bmaps.types()
```
`## [1] "Aerial"           "AerialWithLabels" "Road"`
```bmaps.plot(altalake, type="AerialWithLabels")
```

The next thing we’ll adjust is the zoom level. The zoom level (level of detail) is calculated automatically, but it may be that you’re looking for higher (or lower) resolution. To specify a resolution specifically, use `res=300` (where 300 is the resolution in dpi; useful when exporting figures), or `zoomin=1`, which will use the automatically specified zoom level and zoom in 1 more.

```bmaps.plot(altalake, zoomin=1) # res=300 will also work
```

Step 3: Add overlays

Next we’ll use the `lines()`, `polygon()`, `arrows()`, and `text()` functions we went over earlier to draw on top of the map we’ve just plotted. What’s important here is that we specifically have to tell {rosm} that we don’t want it to project our data (if you plot with `axes=TRUE` you’ll see that `osm.plot()` is not plotting in lat/lon. It’s plotting in EPSG:3857…which you don’t need to understand but may be useful if you do happen to understand it). This sounds intimidating but it’s actually very easy. For larger areas this may lead to your points not lining up, but luckily there are `osm.lines()`, `osm.points()`, and `osm.polygon()` functions that will automatically project your data).

```# plot without projecting
osm.plot(altalake, project=FALSE)
points(altacores\$lon, altacores\$lat, pch=15, cex=0.6)
text(altacores\$lon, altacores\$lat, labels=altacores\$name, adj=c(-0.2, 0.5), cex=0.5)
```

Step 4: Putting it all together

Putting it all together, an example plotting script might like this:

```library(prettymapr)
library(rosm)

# Find bounding box
altalake <- makebbox(50.1232, -122.9574, 50.1035, -123.0042)
# Plot Bing Maps Aerial
bmaps.plot(altalake, res=300, project=FALSE, stoponlargerequest=FALSE)
# Plot overlays
points(altacores\$lon, altacores\$lat, pch=15, cex=0.9)
text(altacores\$lon, altacores\$lat, labels=altacores\$name, adj=c(-0.2, 0.5), cex=0.9)
```

You’ll notice that it still doesn’t have a scale bar or north arrow, and the margins aren’t exactly how we’d like them. It’s possible to eliminate margins manually using `par(mar=c(0,0,0,0), oma=c(0,0,0,0))`, and add a scale bar using the `addscalebar()` function (in the {prettymapr} package), but the easiest way is to use the `prettymap()` function in {prettymapr} to do it all in one step. This might look like the following:

```library(prettymapr)
library(rosm)

altalake <- makebbox(50.1232, -122.9574, 50.1035, -123.0042)
prettymap({bmaps.plot(altalake, res=300, project=FALSE, stoponlargerequest=FALSE)
points(altacores\$lon, altacores\$lat, pch=15, cex=0.9)
text(altacores\$lon, altacores\$lat, labels=altacores\$name, adj=c(-0.2, 0.5), cex=0.9)})
```

There’s tons of options for `prettymap()` that let you customize the north arrow, scale bar etc., which you can find in the {prettymapr} manual.

Using {rcanvec} plot basemaps

```library(prettymapr)
library(rcanvec)
```

Step 1: Find your bounding box

The {rcanvec} package plots based on a bounding box, or an area that you would like to be visible. There’s a few ways to go about doing this, but the easiest way is to visit the Open Street Maps Export page, zoom to your area of interest, and copy/paste the values into `makebbox(northlat, eastlon, southlat, westlon)` from the {prettymapr} package. You can also use `searchbbox("my location name")`, also from the {prettymapr} package, which will query google for an appropriate bounding box. You’ll notice that the bounding box returned by these methods is just a 2×2 matrix, the same as that returned by `bbox()` in the {sp} package.

```altalake <- searchbbox("alta lake, BC")
# or
altalake <- makebbox(50.1232, -122.9574, 50.1035, -123.0042)
```

Back in the days of paper maps, NTS references were used to number maps in an orderly way. For example, the “Wolfville” 1:50,000 scale mapsheet would be referred to as “021H01”, and our mapsheets for Alta Lake in BC would be 092J02 and 092J03. If you don’t know what these are it’s still ok, but this is how the government organizes the data. Take a minute to get familiar with your NTS reference(s), if you feel into it.

```# get a list of NTS references based on a bounding box
nts(bbox=altalake)
```
```## [[1]]
## [1] "092" "J"   "03"
##
## [[2]]
## [1] "092" "J"   "02"```

If you run this command and your bounding box returns more than 4 mapsheets, you’re probably going to want to zoom in. If you’re looking to export your data, you’ll need to have some idea of what these are.

Step 2: Preview your map

{rcanvec} has a method to quickly plot a bounding box: `canvec.qplot()`. We’ll pass our bounding box as the `bbox` argument, and we can use the `layers` argument to customize our layers. If things take too long to plot, you may want to just use `layers="waterbody"` (the “road” layer in particular usually takes a long time to plot). Layers you may want to use are “waterbody”, “forest”, “river”, “contour”, “building” and “road”. Note that the order you put them in will change the appearance of the map.

```canvec.qplot(bbox=altalake)
```

```canvec.qplot(bbox=altalake, layers=c("waterbody", "river"))
```

Step 3: Refine your plotting options (optional)

It’s possible (but not at all necessary) to load layers individually and plot them manually, giving us more control over the appearance of the map. You’ll have to have called `canvec.qplot(bbox=XX)` or `canvec.download(nts(bbox=XX))` before you can load a layer.

```waterbody <- canvec.load(nts(bbox=altalake), layerid="waterbody")
rivers <- canvec.load(nts(bbox=altalake), layerid="river")
forest <- canvec.load(nts(bbox=altalake), layerid="forest")
```

Using the graphic parameters such as `col` and `border`, we can plot our data manually. Note you’ll have to use `xlim` and `ylim` arguments to zoom in. Also, for all calls to `plot()` after the first one, you’ll have to pass `add=TRUE` or it will create a new plot. Because we’re dealing with 2 mapsheets, canvec.load() returns a list of layers, and if you don’t understand that you should pobably stick to using `canvec.qplot()`.

```plot(waterbody[[2]], col="lightblue", border=0, xlim=altalake[1,], ylim=altalake[2,])
```

If you’d still like to use the `canvec.qplot()` function, it’s also possible to build this customization as a “list of lists”, best shown by example:

```plotoptions = list()
plotoptions\$waterbody <- list(col="lightblue", border=0)
plotoptions\$river <- list(col="lightblue")

canvec.qplot(bbox=altalake, layers=c("waterbody", "forest", "river"), options=plotoptions)
```

Step 4: Add overlays

Next we’ll use the `lines()`, `polygon()`, `arrows()`, and `text()` functions we went over earlier to draw on top of the map we’ve just plotted. Unlike {rosm}, {rcanvec} plots in lat/lon natively, so we don’t have to worry about projections.

```canvec.qplot(bbox=altalake, layers=c("waterbody", "forest", "river"), options=plotoptions)
points(altacores\$lon, altacores\$lat, pch=15, cex=0.6)
text(altacores\$lon, altacores\$lat, labels=altacores\$name, adj=c(-0.2, 0.5), cex=0.5)
```

A neat trick is to use the {rosm} package to add a hillshade on top of our map, which we would normally do before plotting our overlays. We’ll have to tell `osm.plot()` not to project its data, since we’re already in lat/lon.

```canvec.qplot(bbox=altalake, layers=c("waterbody", "forest", "river"), options=plotoptions)
```

Step 5: Putting it all together

Putting it all together, an example plotting script might like this:

```library(prettymapr)
library(rcanvec)
library(rosm)

# Find bounding box
altalake <- makebbox(50.1232, -122.9574, 50.1035, -123.0042)
# Plot waterbody, forest, river, road
canvec.qplot(bbox=altalake, layers=c("waterbody", "forest", "river", "road"))
# Plot overlays
points(altacores\$lon, altacores\$lat, pch=15, cex=0.6)
text(altacores\$lon, altacores\$lat, labels=altacores\$name, adj=c(-0.2, 0.5), cex=0.5)
```

You’ll notice that it still doesn’t have a scale bar or north arrow, and the margins aren’t exactly how we’d like them. It’s possible to eliminate margins manually using `par(mar=c(0,0,0,0), oma=c(0,0,0,0))`, and add a scale bar using the `addscalebar()` function (in the {prettymapr} package), but the easiest way is to use the `prettymap()` function in {prettymapr} to do it all in one step. This might look like the following:

```library(prettymapr)
library(rosm)
library(rcanvec)

altalake <- makebbox(50.1232, -122.9574, 50.1035, -123.0042)
prettymap({canvec.qplot(bbox=altalake, layers=c("waterbody", "forest", "river", "road"))
points(altacores\$lon, altacores\$lat, pch=15, cex=0.6)
text(altacores\$lon, altacores\$lat, labels=altacores\$name, adj=c(-0.2, 0.5), cex=0.5)})
```

There’s tons of options for `prettymap()` that let you customize the north arrow, scale bar etc., which you can find in the {prettymapr} manual.

Exporting CanVec data

If you’re like me, you’ll probably want to make a map in some program other than R (QGIS is my personal favourite). You can do this easily using the `canvec.export()` function. The function combines all the NTS sheets you pass into it and exports to a folder of your choosing (with nice human-readable names!).

```canvec.download(nts(bbox=altalake))
canvec.export(nts(bbox=altalake), "canvec_data")
```
```## Writing dsn: canvec_data; layer: residential_area
## Writing dsn: canvec_data; layer: building
## Writing dsn: canvec_data; layer: building_poly
## Writing dsn: canvec_data; layer: tank
## Writing dsn: canvec_data; layer: power_transmission_line
## Writing dsn: canvec_data; layer: contour
## Writing dsn: canvec_data; layer: glacial_debris_undifferentiated
## Writing dsn: canvec_data; layer: elevation_point
## Writing dsn: canvec_data; layer: permanent_snow_and_ice
## Writing dsn: canvec_data; layer: manmade_hydrographic_entity
## Writing dsn: canvec_data; layer: hydrographic_obstacle_entity
## Writing dsn: canvec_data; layer: river
## Writing dsn: canvec_data; layer: waterbody
## Writing dsn: canvec_data; layer: island
## Writing dsn: canvec_data; layer: pit
## Writing dsn: canvec_data; layer: mine
## Writing dsn: canvec_data; layer: domestic_waste
## Writing dsn: canvec_data; layer: lumber_yard
## Writing dsn: canvec_data; layer: aboriginal_lands
## Writing dsn: canvec_data; layer: nts50k_boundary_polygon
## Writing dsn: canvec_data; layer: ski_centre
## Writing dsn: canvec_data; layer: golf_course
## Writing dsn: canvec_data; layer: park
## Writing dsn: canvec_data; layer: trail
## Writing dsn: canvec_data; layer: campground
## Writing dsn: canvec_data; layer: picnic_site
## Writing dsn: canvec_data; layer: wetland
## Writing dsn: canvec_data; layer: named_feature
## Writing dsn: canvec_data; layer: railway
## Writing dsn: canvec_data; layer: runway
## Writing dsn: canvec_data; layer: road
## Writing dsn: canvec_data; layer: junction
## Writing dsn: canvec_data; layer: forest```