Last month, @freakonometrics posted a bunch of articles on his blog about the North pole (in French):

There are some maps on these articles that were produced using ggplot2 excellent package in R. We want to share our code here.

First, let’s get a world map. We’ll use the getMap function from the rworldmap package. I prefer this function to the map one, because it gives more accurate and more up to date borders. I used the same method when I explained how to draw the European Union map.

library(rworldmap) library(dplyr) library(ggplot2) library(geosphere) library(gpclib) # World map worldMap

That code gives the following map:



There is an awesome function that enables us to change the projection for the map in a split second: coord_map, still in ggplot package. Playing with the orientation is easy too.

worldmap

Now, we wanted to highlight regions in Canada, Russia and Greenland lying 3000km or less away from the North Pole. This seems a pretty mere task, but it gets kind of arduous when dealing with multiple polygons representing a single country.

We need a function that gives the coordinates of points on a circle given a radius and an origin. The tricky thing here, is that the North Pole latitude is \(90^\circ\), but there is an infinity of coordinates for the longitude, and using the destPoint function from the geosphere package does not work when providing the couple \((0,90)\). To get around this issue, I cheated a bit, by adding a really small value to the longitude and latitude of the couple \((0,90)\).

# -- # Given the long/lat coordinates of an origin (x) and a radius (radius) in km, # returns the coordinates of 360 points on the circle of center x and radius radius km. # -- # x (numeric vector) : coordinates of the origin of the circle # radius (numeric) : radius of the circle # -- distantCircle

Now, we are going to use the union function from the gpclib package to create a single gpc.poly object for each country we are interested in. It will be easier to deal with such an object. Let's have a look at the coordinates of Canada:

> head(tail(world.df[world.df$region == "Canada",],10),4) long lat group region 2756 -93.72066 77.63433 Canada.29 Canada 2757 -93.84000 77.52000 Canada.29 Canada 2758 -79.26582 62.15867 Canada.30 Canada 2759 -79.65752 61.63308 Canada.30 Canada

As you can see, there are different groups here, because there are 30 polygons defined to draw the borders of Canada in this data.frame. So, for each of these polygons, we are going to extract the coordinates of the intersection with our circle. The following function gives the gpc.poly object for a country:

# -- # From a country name and a data.frame obtained using fortify # returns a polygon of the country # -- # country (string) # map.df (data.frame) : a data.frame obtained with fortify # -- polygonCountry 1){ for(i in 2:length(country.poly)){ country.poly.tot

And here is the function that gives the coordinates of the intersections of the country's polygons and the North Pole circle of radius 3000 km:

# -- # Returns coordinates of the intersection of a country # and a polygon # -- # x (list) : result from polygonCountry() # poly (gpc.poly) : a polygon # -- getIntersection 0){# if there is an intersection # We have to be careful, because there might be multiple polygons poly.intersect

We can apply these functions for Canada, Russia and Greenland. I want the circles to be drawn last so they are not hidden by the polygons. So, I make sure their name in the data.frame begin with a "z". I think there is a more efficient and beautiful way to accomplish this, but I am not aware of it. If you are, please share your method with me in the comment section below.

# Let's bring out parts of Canada, Russia and Greenland 3000km or less away from North Pole countries

For ggplot to be able to fill only the polygons of our three countries, let's add a column in world.df2, where the value is not available for every other country:

world.df2$fill

And finally, we can print the map!

worldmap

Now, if you want to reproduce the animated gif, you can wrap the above code in a function where the argument is the angle:

rotateMap

Now, let's see how to work with gridded data. To illustrate the method, we are going to use some NASA data : GISS Surface Temperature Analysis. There is a link to a txt file for gridded data at the bottom of the page.

First, we need to load the file in the R session.

temp

If we want to have the same colours and the same breaks on the palette, it is not difficult:

breaks

Then, we can add the interval for the surface temperature anomaly:

temp$interval

The last step is to create the plot and wrap the code in a function (to change, once again, the angle):

rotate_map

Before plotting the map, let's analyze this code. We want to plot tiles, so we use the geom_tile() function (that makes sense). Data come from temp, where we can find the x and y coordinates. Since we want the tiles to be filled with a colour corresponding to the surface temperature anomaly, we make sure to use the fill argument. And finally, in order to get the same colours as those on the NASA map, we have to change them manually, using the function scale_fill_manual.

rotate_map(1)

Once again, if you want to create an animated version, it is easy, but this time, it takes quite a long time...

library(animation) saveGIF({ ani.options(nmax = 360) for(i in seq(0,360)){ print(rotate_map(i)) print(i) } }, interval = 0.1, outdir="/Users/ewengallic/Documents/Projets/carte_pole", movie.name = "temp_anomalies.gif")

The result is a huge file (~28 Mb) which can be seen here (I left the legend title in French).