Motivation

A major part of research is based on spatial analyses all over the planet. R delivers a huge toolset with professional tools, especially in statistic. Beside the professional statistic tools, there are only a few ways to plot large-scale spatial data. The basic functions of levelplot() or plot.raster() are simple and good tools for first testing, but they are not developed enough to meet enhanced cartographic requirements. Therefore, the aim of this publication is to show a simple and fast way of plotting spatial data with a standard cartographic projection.

Fundamentals

In this study, the van der Grinten world map projection is introduced and used for the following plots. It is by far the most and widely used projection for large-scale maps. In general this kind of projection is based on a compromise between the Mollweide (equal area) and Mercator (conformal) projection. The projection is mostly used for global or hemispheric maps.

Coding

For this tutorial, we use the function of R package ggplot and a simple world map as background layer to proof that our function works.

library(ggplot2)
library(maptools)
data(wrld_simpl)
map<-wrld_simpl

To do spatial maps, ggplot2 has implemented the function coord_map(). (Note: always characterize xlim and ylim)

ggplot()+theme_bw()+
  geom_polygon(data=map, aes(long, lat, group=group), fill="grey",color="transparent", alpha=0.65)+
   ggtitle("Sample sites")+
   coord_map(projection="vandergrinten",xlim=c(-180,180),ylim=c(-90,90))
## Regions defined for each Polygons

From the plot above, it is visible that the standard axis labels are not working well with the different projections. This is especially visible, if we focus our map on a small extend, for example Europe.

ggplot()+theme_bw()+
  geom_polygon(data=map, aes(long, lat, group=group),                fill="grey",color="transparent", alpha=0.65)+
   ggtitle("Sample sites")+
   coord_map(projection="vandergrinten",xlim=c(-10,50),ylim=c(30,73))
## Regions defined for each Polygons

Since we know from the plots before, we have to create our own axes labels. This is especially a big problem, if you work with data from the pole. I suggest, that we remove all axes and the corresponding labels.

ggplot() + 
               geom_polygon(data=map, aes(long, lat, group = group), fill="transparent",color="gray70")+
               ggtitle(paste("Sample site"))+theme_bw()+scale_color_manual(labels = c("Eoi400", "E280"), values = c("#0073C299", "#A7303099"))+
         annotate("rect", xmin = -181, xmax = 181, ymin = -3, ymax = -15,fill="white")+
         annotate("rect", xmin = -181, xmax = 181, ymin = -0.25, ymax = -30,fill="white")+
         annotate("segment", x = -180, xend = 180, y =80.8, yend = 80.8,
                  colour = "gray40",size = 0.3)+
        annotate("text",x = 187, y = 59,label = paste0(60, "°N"),size=3.5,angle = -54) +
         annotate("text",x = 186, y = 30,label = paste0(30 ,"°N"),size=3.5,angle = -71) +
         annotate("text",x = -90,y = -6,label = "90°W",size=3.5) +
         annotate("text",x = 90,y = -6,label = "90°E",size=3.5) +
        
         annotate("text",x = -174,y = -3,label = "180°W",size=3.5) +
         annotate("text",x = 0, y = -6,label = "0°",size=3.5) +
         annotate("text",x = 176,y = -3,label = "180°E",size=3.5) +
       
         coord_map(projection="vandergrinten",xlim=c(-185,185),ylim=c(-8,80))+
        theme(panel.border = element_blank(),axis.ticks.y = element_blank(),plot.title = element_text(hjust = 0.5),
               axis.ticks.x = element_blank(),
               axis.text.x = element_blank(),## <- this line
               axis.text.y = element_blank(),
               axis.title=element_text(size=0),panel.grid.major = element_line(colour = "gray40",size = 0.3),
               legend.key.size = unit(1, "cm"),panel.background=element_blank())+
         scale_y_continuous(expand=c(1,1),breaks = seq(0, 90, 30))+
         scale_x_continuous(expand=c(1,1),breaks = seq(-180, 180, 90))
## Regions defined for each Polygons

Let’s try the entire globe. Note, based on the underlying shape file, the poles are not plotted because the quality of these features are to low and would destroy the map.

ggplot() + 
               geom_polygon(data=map, aes(long, lat, group = group), fill="transparent",color="gray70")+
               theme_bw()+
        annotate("text",x = 188, y = 45,label = paste0(45, "°N"),size=3.5,angle = -63) +
  annotate("text",x = 184, y = 0,label = paste0(0),size=3.5,angle = 0) +
         annotate("text",x = 188, y = -45,label = paste0(45,"°S"),size=3.5,angle = 245) +
    
         coord_map(projection="vandergrinten",xlim=c(-185,185),ylim=c(-85,85))+
        theme(panel.border = element_blank(),axis.ticks.y = element_blank(),plot.title = element_text(hjust = 0.5),
               axis.ticks.x = element_blank(),
               axis.text.x = element_blank(),## <- this line
               axis.text.y = element_blank(),
               axis.title=element_text(size=0),panel.grid.major = element_line(colour = "gray40",size = 0.3),
               legend.key.size = unit(1, "cm"),panel.background=element_blank())+
         scale_y_continuous(expand=c(1,1),breaks = seq(-90, 90, 45))+
         scale_x_continuous(expand=c(1,1),breaks = seq(-180, 180, 90))
## Regions defined for each Polygons

Conclusion

To sum up, the presented code parts are a first step for cartographic illustrations in R. It is obvious that the huge amount of manual edits indicate that cartographic illustrations with R are not simple, especially for beginners.


Daniel Balting – CC 3.0