GIS in R using the Google Maps API
5 December, 2019
Everyone loves maps. If you're tired of seeing maps without New Zealand in them, why not learn how to build your own? We'll be querying the Google Maps API to build a base map, and then using some publicly available data to visualise the locations and magnitudes of earthquakes measured around Fiji.
1. The data
We'll be using data from the R datasets library. The library loads by default, so you should be able to access the quakes dataset without explicitly bringing it into the global environment. But let's bring it in anyway and have a look:
library(tidyverse) quakes <- quakesglimpse(quakes)
So we have 1000 observations of 5 variables, the variables being latitude, longitude, magnitude of the quake, depth of the quake, and the number of stations that recorded it. Seems pretty simple!
2. Generating the base map
Firstly, let's activate the API with our Google account. After logging in and clicking the "Get started" button here, we then select the "Maps" product. We'll create a project titled "GIS tutorial". If this is your first time using the Google Maps API, you will need to provide your billing information. However, you'll currently also get $200 USD free credit every month; and since the first 100,000 Static Maps queries will cost you $2, you probably won't pay a cent.
Copy the API key that has been generated for you. Let's load the packages and register your Google API key. register_google() is a function from the ggmap package. Finally, remember to keep your API key private!
library(tidyverse)
library(ggmap)
register_google(key = "[your_key_here]")
Next, we'll create our base map with ggmap::get_googlemap(). There are a few important parameters here.
We pass longitude and latitude coordinates to the `center` parameter, which is to be the center of the map. To get those coordinates, we'll just do a simple approximation of somewhere in the middle of our data:
c(median(quakes$long), median(quakes$lat))
[1] 181.41 -20.30
So we have a central longitude of 181.41. Trying to pass this to get_googlemap won't work, however, because they're plotting points with a coordinate system that has longitudes between -180 and 180. Our longitude of 181 is their longitude of -179! So let's add a new column to our quakes tibble, containing a representation of longitude that Google Maps will be happy with. If the current longitude value is over 180, we convert it to a negative representation. Our new longitude column is called... new_long! Of course, our new central longitude of 181.41 becomes -360 + 181.41 = -178.59. Note: For some reason, Google Maps is not happy with a mixture of negative and positive longitude values crossing the meridian, but it's fine with all negative values... so for this reason, we'll just convert everything to a negative longitude:
quakes <- quakes %>%
mutate(new_long = -360 + long)
Moving on to the rest of the get_googlemap() parameters! We can adjust the zoom using the `zoom` parameter. (Unfortunately, only integer values are supported, which means that we may need to use another way of obtaining the most appropriate zoom for our map.) We can also use the `maptype` parameter to select from terrain, satellite, road and hybrid map types. The default is terrain, which is what we're after here, so we don't need to explicitly specify that. If we feel like it, we can specify a black-and-white output by passing "bw" to the `color` parameter, but let's just leave it at the default colour now.
Let's create our `ggmap` map object, and then plot it.
map <- get_googlemap(
center = c(-178.59, -20.30),
zoom = 6)
ggmap(map)
Great! But what if we want more control over the layers of objects that are plotted over the base map? We can pass
a string containing a bunch of style settings to the `style` parameter. The syntax for this is rather
obtuse. Each `style` declaration can optionally contain a `feature` (e.g., a road),
an element (that is, a characteristic of a feature, e.g., a road's name), and, finally, the required
style rules. I encourage you to read more about it in the official documentation. In this case, we're going to remove road markings, labels and administrative markings from our map.
Referring back to the style nomenclature, to remove a feature (road), we set its element (visibility) "off".
map <- get_googlemap(
center = c(-178.59, -20.30),
zoom = 6,
style = "feature:road|visibility:off&style=element:labels|visibility:off&style=feature:administrative|visibility:off")
ggmap(map)
3. Plotting
Now that we've built our base map, we don't have to call the Google API again unless we need to change the map. We can just apply normal ggplot geoms over the ggmap() call to our map object! Let's do that, by plotting the earthquake locations with geom_point, and colouring by magnitude.
ggmap(map) +
geom_point(data = quakes, aes(x = new_long, y = lat, colour = mag))
See that warning message, telling us it's removed 389 rows containing missing values? That's because our map is too zoomed-in to display all the points in the dataset. We actually have to set the zoom parameter to 4, in order for us to plot every point on the map. But that's really zoomed out! Fortunately, we can specify the limits of the x- and y-axes using `coord_fixed()`, and this acts as a zoom function that we can fine-tune.
To do this, we pass longitude and latitude limits to `xlim` and `ylim` parameters, respectively. We also need to specify the aspect ratio so that our map is not warped. We pass our latitude value into the `ratio` parameter of `coord_fixed()`. Because our latitude = -20.30, we pass the following:
ratio=1/cos(pi*--20.30/180))`.
For more information, see user Pere's excellent answer in this Stack Overflow post.
Let's see this in action. Because we've already created our base map above (and stored it as `map`), we don't need to repeat ourselves. We simply plot the map, and apply a new coordinate scheme. I've chucked in some lat/lon limits just by having a look at the map and seeing what should work.
ggmap(map) +
coord_fixed(xlim = c(-170, -200), ylim=c(-40, -10), ratio=1/cos(pi*--20.30/180)) +
geom_point(data = quakes, aes(x = new_long, y = lat, colour = mag))
To be continued...