R provides a rich environment for cartographic plotting.

---

Maps with ggplot2

The powerful ggplot2 library is capable of producing quite nice looking maps. In order to showcase these capacities we will visualize the location of weather stations maintained by Deutscher Wetterdienst (DWD) (German Weather Service). The data was downloaded from the DWD data portal on April 21, 2017.

First, before plotting any data we have to prepare our data set. Hence, we start by importing the data set and assign it a proper name. In this example we use the read.csv2() function for loading the data.

dwd <- read.csv2("https://userpage.fu-berlin.de/soga/data/raw-data/DWD.csv",
  stringsAsFactors = FALSE, sep = ","
)

The dwd data set consists of 599 rows, each of them representing a particular weather station in Germany, and 22 columns, each of them corresponding to a variable/feature related to that particular weather station. For the purpose of this tutorial we are only interested in the variables LAT and LON, corresponding to the geographic coordinates of each particular weather station in the data set, as well as the record length in years (RECORD.LENGTH). Hence, we subset our data set based on these variables. Further, we make sure that we exclude all missing values.

## subset data
dwd <- dwd[, c("LAT", "LON", "RECORD_LENGTH")]
## drop missing values
dwd <- dwd[complete.cases(dwd), ]

dwd$LAT <- as.numeric(dwd$LAT)
dwd$LON <- as.numeric(dwd$LON)
nrow(dwd)

## [1] 599

head(dwd)

##       LAT     LON RECORD_LENGTH
## 1 47.8413  8.8493            55
## 2 50.7827  6.0941           160
## 3 52.9335  8.2370            45
## 4 48.2156  8.9784            30
## 5 48.6159 13.0506            64
## 6 52.4853  7.9126            55

For the purpose of visualization we categorize the variable RECORD.LENGTH into five groups:

• group 1: < 10 years
• group 2: 10-30 years
• group 3: 30-60 years
• group 4: 60-90 years
• group 5: > 90 years

To achieve this task we apply the mutate() function provided by the dplyr package. If you get an error message, type install.packages('dplyr') to install the package.

library(dplyr)
dwd <- dwd %>%
  mutate(RECORD_LENGTH_CATEGORY = cut(RECORD_LENGTH,
    breaks = c(-Inf, 10, 30, 60, 90, Inf),
    labels = c(
      "very short (<10)",
      "short (10-30)",
      "middle (30-60)",
      "long (60-90)",
      "very long (>90)"
    )
  ))

Now we are ready to plot a map to visualize the spatial distribution of the DWD weather stations. For this we rely on the raster package and on the ggplot2 package.

library(raster)
library(ggplot2)

If you get an error message saying that one or both of the packages are not found, type install.packages(c('raster', 'ggplot2')) to install the packages.

Basically, we only use the getData() function from the raster package. This function provides us with geographic data for Germany. Type help(getData) into your R console for further details.

# Retrieve Federal States by the the getData() function from the raster package
germany <- getData(country = "Germany", level = 1)

## Warning in getData(country = "Germany", level = 1): getData will be removed in a future version of raster
## . Please use the geodata package instead

## Please note that rgdal will be retired during October 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-7, (SVN revision 1203)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.6.2, released 2023/01/02
## Path to GDAL shared files: C:/Users/mceck/AppData/Local/R/cache/R/renv/cache/v5/R-4.3/x86_64-w64-mingw32/rgdal/1.6-7/10b777236c9e7855bc9dea8e347e30b7/rgdal/gdal
##  GDAL does not use iconv for recoding strings.
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 9.2.0, March 1st, 2023, [PJ_VERSION: 920]
## Path to PROJ shared files: C:/Users/mceck/AppData/Local/R/cache/R/renv/cache/v5/R-4.3/x86_64-w64-mingw32/rgdal/1.6-7/10b777236c9e7855bc9dea8e347e30b7/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-1
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.

In the next step we build a ggplot2() graph, layer by layer.

• Firstly, we initialize the plot (ggplot()).
  
• Then we add a geom layer, more specifically a geom_polygon(). Here, we specify the variables to find longitude and latitude as well as defining the color of the borders and the filling of the polygon.
  
• Then we add another geom layer (geom_point()), representing the locations of the DWD weather stations. Again, we specify latitude and longitude, as well as color and size parameters.
• Then we provide the coordinate system (coord_map()).
  
• Finally, we set some layout aesthetics (theme_bw(),xlab(),ylab(),labs() and ggtitle()).

library(mapproj)

ggplot() +
  geom_polygon(
    data = germany,
    aes(x = long, y = lat, group = group),
    colour = "grey10", fill = "#fff7bc"
  ) +
  geom_point(
    data = dwd,
    aes(x = LON, y = LAT, col = RECORD_LENGTH_CATEGORY),
    alpha = .5,
    size = 1.5
  ) + 
  coord_map() +
  theme_bw() +
  xlab("Longitude") +
  ylab("Latitude") +
  labs(color = "Record length") +
  ggtitle("DWD Weather Stations")