Maps & geospatial data

Everyone loves a map. Whether you’re publishing a blog post, an academic article, or a data visualization, they’re (almost) always a great way to show off your skills. However, they can be a major pain to make and get right.

Let’s get right into it!

Getting country boundaries the easy way with `rnaturalearth`

A nice, easy way of getting country or region boundaries is the rnaturalearth package. Natural Earth is a great resource for this kind of data, and this package is just an easy way to import it.

One function imported with this package is the ne_countries() function, which allows you to get country boundaries. You can specify the scale of the data you want, and the continent you want to get. If you don’t specify a continent, you’ll get the whole world.

hHe data comes in different levels of detail, and the 110 scale is the least-detailed version you can get, but it’s good enough for our use case. (There is also 50 and 10 scale)

If you want high resolution map data, you’ll have to download the rnaturalearthhires package from Github.

You can do this with the devtools package, which installs the package from Github.

But what are we looking at here? This is just a normal data frame that includes one special column, ‘geometry’. Each geometry value is a polygon depicting the shape of the country.

We also, as a nice bonus, get a whole bunch of other data from Natural Earth, such as the population, name in many different languages, ISO codes, and region.

Generally, you don’t want to touch these geometry columns, just let them do their own thing behind the scenes. They’re just big lists of numbers that describe the shape of the country.

Making maps with `sf()`

So how do we plot this geometry to a map?

With the sf package loaded, your you have access to some new functions, including here the geom_sf() function, which allows you to plot the data in ggplot. You don’t need to specify the aes() function to show basic shapes.

Just like that, we’ve made a map!

Working with cartesian coordinates

You’ll notice that the X and Y coordinates are a little different than what we’re used to: they are in decimal degrees. Whereas a lot of real-world maps work in degrees, minutes, and seconds, with 60 minutes in a degree, and 60 seconds in a minute, these maps are in the base-10 equivalent.

For example, New York City is at 40°42′46″North, 74°0′22″ West, or 40.712778, -74.006111 in decimal degrees. Google Earth is a good place to find these, the longitude and latitude will show up in the lower right corner.

In the settings, you might have to change it to decimal degrees.

If you need the coordinates of a location, Wikipedia is also a pretty good resource.

It’s easy to get mixed up about Longitude and Latitude and X and Y, so here are a couple tricks:

Remember that Lats lay Flat!

And that Y has 3 lines, lat has 3 letters!

And that Long has the long lines! (and the larger numbers)

All longitude lines are the same length, and they all meet at the poles. Latitude lines are shorter as you get closer to the poles.

So to crop our map, we can simply set the limits of the x and y axes with lims(), just like in any other ggplot object. Let’s focus on just Southeast Asia for this segment. Because we’re still using ggplot, which thinks in terms of X and Y, we’ll have to use those to set the limits.

Classwork: Crop a map

1. Choose a country.
2. If it’s not in Asia, download the data for it using ne_countries().
3. Manually crop the map to that country and its neighbors, but loosely so that you can see the countries around it.
4. Make the map square-ish.

For example, I made this nice map of Myanmar:

Layering data from different objects

So far, we’ve been largely making our plots from one data frame at a time, but in some cases, we’ll want to layer things on top of each other. For example, what if I want to add the different states and regions of Myanmar (it has both!), and just use this as a background map?

I could first use the ne_states() function to get the states of Myanmar:

Alternatively, if you don’t have the naturalearthhires package installed, you can use the ne_download() function to download the data directly from Natural Earth. You can then just filter for Myanmar.

I would then add both to a blank ggplot() object, and then add the two layers on top of each other. I would indicate which data set to use for each layer using the data= argument. Each layer can have its own aes() function, so we can style each one of them separately.

Labeling a map the easy way with `geom_sf_label()`

Nice! But maybe our reader’s geography is terrible, and they don’t know their Ayeyarwadys from their Tanintharyis, so we might have to add some labels to the map. The built-in way to do this is with geom_sf_label() or geom_sf_text(), which is a lot like geom_label() but works with the geometry column.

Your first reaction might be “wow, that’s some terrible overlapping, and some countries are missing!” and you’d be partially correct.

Remember that the map you see on screen is just a little preview of what gets exported through ggsave(), and you should always look at the final project. By saving it and opening it in an image editor, you can see the final result. Do this before you make any irrational decisions.

If we have time at the end of class, I’ll teach you a better, but more complicated, way to label a map.

Classwork: Make a map

1. Add your country of interest to the map you made previously.
2. Make it a separate color from the background map, and fade the colors in the background.
3. Use the following code to download the capitals of every country:

4. Add the capitals to your map as as both points and labels.
5. Move the labels around using vjust and hjust, so they don’t overlap with the points. For example, this would would vertically

6. Save it as a PNG file, and review the final project

Your result might look something like this:

Cartography can get you in a lot of trouble

Wherever you publish your data, you have to be very careful about the maps you present; some countries are very touchy about this.

This can apply to your own country, too, and can get you in some deep trouble. Wherever you publish, you’ll have to consider:

1. The laws of the country you’re in,
2. The laws of the country you’re publishing in,
3. The laws of the country you’re from,
4. The laws of any country you might want a visa for in the future,
5. Your own politics and viewpoints,
6. The politics and viewpoints of the people you’re publishing for,
7. How many angry comments you want to get on social media,
8. etc, etc, etc.

This might be completely untintentional: In the above map of Myanmar, the border between India and China is disputed, and we’ve just inadvertently taken a position on it, making this map illegal to publish under Chinese law.

This can be a real headache, and if this worries you, you might not want to publish maps at all.

Getting country data from shapefiles

These download functions are nice, but maybe you’ll need to need to get data from some other source. The most common format for these are shapefiles, which are a collection of files that describe the shape of a country. These are often used in specialized GIS software, but can be imported into R using the sf package.

For example, if you want to avoid the can of worms described above, naturalearthdata.com has shapefiles from each country’s point of view available to download directly.

The default shapefile from ne_countries() depicts Crimea as belonging to Russia, but let’s say you want to show a map depicting it as Ukrainian. You’ll want to:

1. Download the .zip file, from, say, the Swedish perspective.
2. Unzip it.

You’ll notice that you get 6 files, not just one. What a deal!

The other folders are just metadata and database files that are necessary for the shapefile to work, and you can ignore them for now. The only one we care about is the .shp file, which contains the actual data.

But if you move the file, be sure to move the other ones too, or the whole thing will break!

Now, when you make a map and publish it somewhere, you’ll have an entirely different group of people angry at you in the comments sections!

Adding data from point coordinates

Sometimes, our data isn’t in a shapefile, but in a CSV or some other format. In this case, we need to do a few things to get it into a format we can use. I found here a nice .csv of solar power plants in France, from here which includes the coordinates of each plant.

To put this on our map properly, we need to convert the longitude and latitude columns into a geometry column. We can do this using the st_as_sf() function, which takes a data frame and converts it into a spatial object.

We also need to set the coordinate reference system (CRS), which is the system that describes how the coordinates are defined. The most common one is WGS84, so we’re just going to hope it’s correct.

Now we can add the power plants to our map using the geom_sf() function, just like we did with the countries.

I’ve also cropped it to just metropolitan France, so we can see the power plants better.

Ok, not a bad start, but maybe a little bit too busy.

Because there are so many solar power plants, maybe we can handle them a little differently, and make a cloropleth map of the number of plants in each region. We can do this by grouping the data by region, and then counting the number of plants in each french department.

Spatial computation

Preparing our data

First, we need to get the department boundaries of France, which I’ll once again take from Natural Earth.

We only want some administrative codes, the names, and the geometry, so I remove everything else.

Next, I’m going to make a temporary data frame with the geometry and the row id of the solar power plants.

Some of these data sets get absolutely massive, so I try to keep my data frames as slim as possible.

Perfoming a spatial join with `st_join()`

The sf library also comes with a variety of functions starting with “st_”, which are used to perform spatial operations. I can’t find it in the documentation anywhere, but I assume it stands for “spatial transformation.” One of these is st_join(), which allows you to join two data frames based on their location.

Just like a normal join, you can use it to combine two data frames based on a common feature. In this case, we want to join the solar power plants with the department boundaries, so we can count the number of solar power plants in each department.

In this case, we use st_join to tell us which department of France each solar power plant is in.

What we get is a data frame with the geometry of the department, and the row id of the solar power plant.

Dropping our geometry layer with `st_drop_geometry()`

We no longer need the geometry of the solar power plants, so we can drop it using the st_drop_geometry() function. This will remove the geometry column from the data frame, and leave us with just the department boundaries and the row id of the solar power plants.

It is still a spatial object, so we can turn it back into a normal data frame with as_tibble().

Now we can simply use the count() function to count the number of solar power plants in each department.

We want to use this data to color our cloropleth map, so we need to join it back to the department boundaries. We can do this using the normal old left_join() function, which will add the number of solar power plants to each department.

I’ve also renamed the column to something a little bit more descriptive.

We can see that this still leaves us with some NA values, which we can clean up using replace_na()

When we plot this, we can see that we’ve successfully counted up the number of solar power plants in each department.

There’s a lot, lot more to talk about with this sort of things, and the whole field of GISc is devoted to it. Things get really complicated really fast, and R isn’t necessarily the best tool to solve the crazy complicated problems. However, it is good enough for 90% of use cases.

Combining plots

But somebody in Saint Pierre and Miquelon is already shouting “Tabarnouche! We are part of France too!” and they are right. So let’s add the overseas territories to the map. We can do this with the cowplot library, which allows you to combine multiple ggplot objects into one.

We need to make separate plots for some overseas territories, and save them as different environment variables. But how can we do this?

I would first save a map with no limits set, and also the legend and title removed. We do this to have as much of the map as possible in one block of code, so that if we have to, say, change the colors, we don’t have to edit every single sub-map.

I’ve also added a scale bar in the top right using the annotation_scale(location = 'tr') function from the ggspatial package.

I would then make the main map, and add the legend, title, etc. to that one.

You’ll get a warning saying “Scale on map varies by more than 10%, scale bar may be inaccurate,” and this is correct.

This is because we’re using the Mercator projection, which distorts the size of the countries as you get further from the equator. However, this is still good to include to give viewers a rough ballpark of how big the different territories are.

I would then make separate plots for each overseas territory I wanted to include.

For the sake of this example, I will just choose 4. I make them all the same dimensions (no super long or tall maps), so that the final map doesn’t look so wonky.

Next, we can combine the plots using the plot_grid() function from the cowplot package. This allows us to arrange the plots in a grid, and specify how many columns we want. For example, I want all the overseas territories in one column, so I set ncol = 1.

Finally, we add this to the side of the main plot. Because we want the main plot to be larger, we fiddle around with the relative widths of the two plots using rel_widths. Each plot is assigned a width, with a larger number meaning a wider plot.

Again, because we only care about the final output, we can save this to a file using ggsave(), and only examine the final output.

Extra: coordinate systems

By default, sf uses the WGS84 coordinate system, which plots the data onto a rectangular grid. However, while this is very useful for navigating a boat across the ocean, it is not always the best way to plot data on a map; locations near the poles are too large and too flat, and distances are heavily distorted.

Using the coord_sf() function, you can set the coordinate system to something more appropriate for your data. For example, this code re-projects our world map to the Robinson projection, which abandons the rigid grid of the Mercator and better preserves the shape of the continents.

Or here’s the UN map (a Polar azimuthal equidistant projection).

All these are done with special codes passed to the coord_sf() function. You can find a list of all the different projections here.

Some countries have an official projection, which might be nice to use for your maps. For example, Switzerland has theSwiss Oblique Mercator.

A useful tool for generating these codes can be found at projectionwizard.org.

Extra: An ultra-handy debugging tip

The largest problem with cartography is the issue of projecting spherical coordinates onto a flat plane. R can sometimes struggle with this, ans has two different methods of doing so. If you’re getting weird errors, try switching between the two methods using these commands:

You can read about what they mean yourself if you want to.

Extra: Map labeling solution No. 1

Back to our map of Myanmar!

Labeling our map correctly is honestly really hard sometimes. The easiest, and often the best way, is simply to export the map as an .svg file, and edit it in a vector graphics editor like Inkscape (free!) or Adobe Illustrator.

While I’m generally an advocate for the “spend 2 hours automating a task rather than spend 1 hour doing it manually” approach to computing, this is one of those cases where it really is easier to just do it by hand.

Extra: Limits or cropping

If you want to crop the map to a specific area, you can use the st_crop() function. This will crop the map to a specific bounding box, which is defined by the minimum and maximum x and y coordinates.

This gets rid of all the data not within that bounding box, which can be helpful in many situations.

Extra: Map labeling solution No. 2

Ok, you say, but I really want to do this in R. I get it.

There exists a great library for labeling called ggrepel. Using its geom_label_repel() function, we can add labels to anything on a plot. If there are too many labels in an area, it will automatically move them around so they don’t overlap. If there are still too many, it will leave some unlabeled.

However, it doesn’t work with geom_sf_label(), so we have to do a little bit of cheating.

First, we need to get the center points (centroids) of each country in the asia data frame. We can do this using the st_centroid() function, which will give us the center of each country.

We use st_centroid() on the asia_cropped data frame, which is the cropped version of the asia data frame. This will give us the center of each country in the cropped area. This prevents our centroid for China and India from being in the center of the country, but instead in the center of the part of that country seen on the map. This will keep the labels on our map.

Next, we need to get the coordinates of each centroid, which we can do using the st_coordinates() function.

st_coordinates() returns a list of X Y, (and maybe Z: elevation and M: measure), so we need to use some janky code to get the longitude and latitude for each centroid.

Now, assuming we haven’t touched our coordinate system at all, we can overlay the centroids on top of the map using geom_label_repel().

Practice and Homework

1. Make a map of your choice. It can be as simple or as complicated as you want.
2. Please email me the code you used to make your map in a document named week_9_homework_(your_name).R by Tuesday, April 6th.
3. Please also include your final map as a PNG or JPG file.