updated world mapping vignette

Former-commit-id: d0b8e9f

vignettes/world-great-circles.Rmd

@@ -1,49 +1,118 @@
-Great circles lines on a world map
+Great circles lines on a world map with **rworldmap** and **ggplot2**
 ========================================================
 
 Sometimes you will want to plot maps on a much larger
-scale that we have covered in the majority of
-the [Creating-maps-in-R repository](https://github.com/Robinlovelace/Creating-maps-in-R/).
+scale that we have covered previously in the
+'Introduction to visualising spatial data in R'
+[tutorial](https://github.com/Robinlovelace/Creating-maps-in-R/raw/master/intro-spatial-rl.pdf),
+hosted on the
+[Creating-maps-in-R github repository](https://github.com/Robinlovelace/Creating-maps-in-R/).
 For this there are a range of options, including packages called
-[*maps*](http://cran.r-project.org/web/packages/maps/index.html),
-[`map_data`](http://docs.ggplot2.org/0.9.3/map_data.html) from the
-*ggplot2* package and [*rworldmap*](http://cran.r-project.org/web/packages/rworldmap/index.html).
-
-We will use the latter two options to show how maps of the entire world
-can easily be produced in R. Amazingly, in each of the packages, the geographic
-data for the world and many of its subregions have already been saved, saving the
+[**maps**](http://cran.r-project.org/web/packages/maps/index.html),
+a function called [`map_data`](http://docs.ggplot2.org/0.9.3/map_data.html) from
+**ggplot2** package and [**rworldmap**](http://cran.r-project.org/web/packages/rworldmap/index.html).
+
+In this post we will use the latter two (newer) options
+to show how maps of the entire world
+can easily be produced in R and overlaid with shortest-line
+paths called *great circles*. Amazingly, in each package, the geographic
+data for the world and many of its subregions
+are included, saving the
 need to download and store files of unknown quality from the internet.
 
-# Plotting continents and great circle lines in base graphics
+## Plotting continents and great circle lines in base graphics
 
+The first stage is to load the packages we'll be using:
 
 ```{r}
-library(rworldmap)
-library(geosphere)
+x <- c("rworldmap", "geosphere", "ggmap")
+lapply(x, require, character.only = T)
+```
 
-s <- getMap()
-plot(s)
-bbox(s)
+Let us proceed by loading an entire map of the world from
+the **rworldmap** function `getMap`:
+
+```{r, fig.keep='none'}
+s <- getMap() # load the map data
+class(s) # what type of are we dealing with?
+nrow(s) # n. polygons
+plot(s) # the data plotted (not shown)
+bbox(s) # the bounding box... of the entire world
+```
+
+The above shows that in single line of code we have loaded
+`s`, which represents the entire world and all its countries.
+This impressive in itself,
+and we can easily add further details like colour based on
+the countries' attributes (incidentally, you can see
+the attribute data by typing `s@data`).
+
+## Adding points randomly scattered over the face of the Earth
+
+But what if we want to add up points to the map of
+the world and join them up? This can be done in
+the same way as we'd add points to any R graphic.
+Using our knowledge of `bbox` we can define the limits
+of random numbers (from `runif`) to scatter points randomly
+over the surface of the earth in terms of longitude. Note the use of
+`cos(abs(l))` to avoid oversampling at the poles,
+which have a much lower surface area than the equator, per
+[line of longitude](http://en.wikipedia.org/wiki/Cylindrical_equal-area_projection).
 
+```{r}
 set.seed(1984)
 n = 20
 x <- runif(n=n, min=bbox(s)[1,1], max = bbox(s)[1,2] )
-y <- runif(n=n, min=bbox(s)[2,1], max = bbox(s)[2,2] )
+l <- seq(from = -90, to = 90, by = 0.01)
+y <- sample(l, size = n, prob = cos(abs(l) * pi / 180))
 p <- SpatialPoints(matrix(cbind(x,y), ncol=2), proj4string=CRS("+proj=longlat +datum=WGS84"))
 
+plot(s)
 points(p, col = "red")
+```
+
+## Joining the dots
+
+So how to join these randomly scattered points on the planet?
+A first approximation would be to join them with straight lines.
+Let's join point 1, for example, to all others to test this method:
+
+```{r}
+plot(s)
+segments(x0 = rep(coordinates(p[1,])[1], n), y0 = rep(coordinates(p[1,])[2], n),
+         x1 = coordinates(p)[,1], y1 = coordinates(p)[,2])
+```
 
-head(gcIntermediate(p[1,], p[2])) # take a look at the output of the gcIntermediate function
-lines(gcIntermediate(p[1,], p[2,]))
+(Incidentally, isn't the use of `segments` here rather clunky - any suggestions
+of a more elegant way to do this welcome.)
+The lines certainly do join up, but something doesn't seem right in the map, right?
+Well the fact that you have perfectly straight lines in the image means bendy
+lines over the Earth's surface: these are not the shortest,
+[great circle](http://en.wikipedia.org/wiki/Great_circle) lines.
+To add these great circle lines, we must use the **geosphere** package:
+
+```{r}
+head(gcIntermediate(p[1,], p[2]), 2) # take a look at the output of the gcIntermediate function
+plot(s)
+lines(gcIntermediate(p[1,], p[2,]), col = "blue", lwd = 3)
 
 # for loop to plot all lines going to zone 5
 for(i in 1:length(p)){
-  lines(gcIntermediate(p[5,], p[i,]), col = "green")
+  lines(gcIntermediate(p[1,], p[i,]), col = "green")
 }
 ```
 
+Fantastic. Now we have great circle lines represented on a
+map with a [geographic coordinate system (CRS)](http://en.wikipedia.org/wiki/Geographic_coordinate_system)
+(as opposed to a projected CRS, which approximates Euclidean distance).
+
+## Beautifying the map
+
+The maps we created so far are not exactly beautiful.
+Let's try to make the map look a little nicer:
+
+
 ```{r}
-# beautify the map a little
 names(s@data)
 library(rgdal)
 # s <- spTransform(s, CRSobj=CRS("+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"))
@@ -58,7 +127,7 @@ for(i in 1:length(p)){
 par(bg = 'white')
 ```
 
-# Doing it in ggplot2
+## Doing it in ggplot2
 
 The 'beautified' map above certainly is more interesting visually, with added
 colours. But it's difficult to call it truly beautiful. For that, as with
@@ -72,8 +141,9 @@ m
 ```
 
 When we add the lines in projected maps (i.e. with a Euclidean coordinate system)
-based solely on origins and destinations, this works fine, but generates incorrect
-shortest path lines in ggplot2.
+based solely on origins and destinations, this works fine, but
+as with the previous example, generates incorrect
+shortest path lines:
 
 ```{r}
 # adding lines
@@ -85,13 +155,15 @@ shortest path lines in ggplot2.
 p1 <- coordinates(p[5,])[rep(1, n),]
 p2 <- coordinates(p)
 
-ggplot() + geom_segment(aes(x = p1[,1], y = p1[,2], xend = p2[,1], yend = p2[,2]))
+# test plotting the lines
+# ggplot() + geom_segment(aes(x = p1[,1], y = p1[,2], xend = p2[,1], yend = p2[,2]))
+
 ggplot() + geom_polygon(data = s,aes(x=long, y=lat, group=group), 
   fill="green", colour="black") +  
   geom_segment(aes(x = p1[,1], y = p1[,2], xend = p2[,1], yend = p2[,2]))
 ```
 
-# Adding great circle lines
+## Adding great circle lines to ggplot2 maps
 
 Adding great circle lines in ggplot2 is similar, but we must
 save all of the coordinates of the paths in advance before plotting,
@@ -125,7 +197,6 @@ empty ggplot2 instances, ready to be filled with layers.
 This is more flexible than stating the data at the outset.
 
 ```{r}
-ggplot() + geom_path(data = paths, aes(lon, lat , group = group))
 ggplot() + geom_polygon(data = s, aes(x=long, y=lat, group=group), 
   fill = "green", colour="black") +
   geom_path(data = paths, aes(lon, lat , group = group)) +
@@ -150,12 +221,28 @@ m + coord_map()
 # remove fill as this clearly causes problems:
 m <- ggplot() + geom_path(data = s, aes(x=long, y=lat, group=group), colour="black") +
   geom_path(data = paths, aes(lon, lat , group = group)) 
-m + coord_map("bicentric", lon = 0)
-m + coord_map("bonne", lat= 0)
+# m + coord_map("bicentric", lon = 0)
+# m + coord_map("bonne", lat= 0)
 m + coord_map("ortho", orientation=c(41, -74, 0)) # for ortho maps
 ```
 
-
+## Conclusion
+
+We've seen 2 ways of plotting maps of the world and overlaying
+'great circles' lines on them. There are probably more, but
+these two options seem to work well, except with
+the bugs in **ggplot2** for plotting polygons in
+many map projections. The two methods are not incompatible
+(see `fortify` for plotting **sp** objects in **ggplot2**)
+and can be combined in many other ways.
+
+For more information on plotting spatial data in R,
+I recommend checking out R's range of
+[spatial packages](http://cran.r-project.org/web/views/Spatial.html).
+For an introductory tutorial on visualising spatial data
+in R, you could do much worse than start with
+[Visualising Spatial Data in R](https://github.com/Robinlovelace/Creating-maps-in-R/raw/master/intro-spatial-rl.pdf)
+by [James Cheshire](http://spatial.ly/) and [myself](http://robinlovelace.net/).