Add a relief (GEOtiff, .tif) to a ggplot of the swiss country borders ((polygon) shapefile, .shp)

Question

I plotted the area of Switzerland (polygon shapefile) and added points (swiss weather stations) via coordinates.

# Boundaries with data-points plotted

library(rgdal)
library(readxl)
library(sp)

library(ggplot2)
library(maptools)
library(plyr)
library(raster)

# import swiss country frontiers (.shp file)

gb <- readOGR("swissBOUNDARIES3D_1_3_TLM_KANTONSGEBIET.shp")

# import coordinates of weather stations from excel file (.xlsx)

coord <- read_excel("SMN-Stationen_20151222.xlsx")
head(coord)
# A tibble: 6 x 10
    SINCE_DT NAT_IND_TX NAT_ABBR_TX           NAME_TX X_KM_COORD_NU     Y_KM_COORD_NU HEIGHT_ASL_NU
      <dttm>      <chr>       <chr>             <chr>         <dbl>             <dbl>         <dbl>
1 2015-12-15       0600         ARO             Arosa        771030            184826          1878
2 2015-10-27       3420         LAC Lachen / Galgenen        707637            226334           468
3 2015-09-08       8040         VEV  Vevey / Corseaux        552106            146847           405
4 2015-08-25       7770         MAR    Les Marécottes        567375            107577           990
5 2015-08-18       5290         FRU          Frutigen        616926            160532           753
6 2015-04-28       0380         BIV             Bivio        771282            148120          1856
# ... with 3 more variables: LONGITUDE_NU <dbl>, LATITUDE_NU <dbl>,     NAME_G_TX <chr>

# plot with ggplot2

gbb <- fortify(gb)

ggplot() +
  geom_polygon(data = gbb, aes(long, lat, group = group, fill =         c("grey40")), color = "white") +
  geom_point(data = coord, aes(x = X_KM_COORD_NU, y = Y_KM_COORD_NU), colour     = "black", size = 3) +
  geom_point(data = coord, aes(x = X_KM_COORD_NU, y = Y_KM_COORD_NU), colour     = "skyblue1", size = 1.5) +  
  scale_fill_identity()

plot image

Now I would like to underlay the graphic with the topographical relief, which I have in form of a GEOtiff file.

When I try to plot the raster the error: "cannot allocate vector of size 196.7 Mb" shows up.

An other problem is, that the two files don't have the same extent. So I can't crop them by the extent. Is there maybe a way to overlay via their CRS?

relief <- raster("25_HYPSO.tif")
> extent(relief)
class       : Extent 
xmin        : 2317050 
xmax        : 3057050 
ymin        : 912950.2 
ymax        : 1412916 
> projection(relief)
[1] "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1     +x_0=2600000 +y_0=1200000 +ellps=bessel +units=m +no_defs"
> 
> extent(gb)
class       : Extent 
xmin        : 485411 
xmax        : 833840.7 
ymin        : 75269.68 
ymax        : 295934 
> projection(gb)
[1] "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1     +x_0=600000 +y_0=200000 +ellps=bessel +units=m +no_defs"

You can download the data from the Federal Office of Topography "swisstopo".

Countryborders (swissBOUNDARIES3D_1_3_TLM_KANTONSGEBIET.shp) can be found in this map:

https://shop.swisstopo.admin.ch/de/products/landscape/boundaries3D

Topographical relief (25_HYPSO.tif) can be found in this map:

https://shop.swisstopo.admin.ch/de/products/maps/national/vector/smv1000

The end product should look something like Timo Grossenbacher did in his blog post (https://timogrossenbacher.ch/2016/12/beautiful-thematic-maps-with-ggplot2-only/). I tried his approach but did not find success because he edited his relief file in advance and I don't know how.

example image, Timo Grossenbacher

score 7 · Accepted Answer · edited Aug 12 '19 at 01:42

I managed to do it now. Timo Grossenbacher told me that he changed the CRS (from CH1903+/LV95 to CH1903/LV03) and clipped (exporting the extent of the swiss country border) with the open source software QGIS.

So that was what I did with the same .tif file but an other .shp file for the export which only contains the country borders (/swissBOUNDARIES3D_1_3_TLM_LANDESGEBIET.shp, can be found in the same download).

I called the resulting file: "landesgebiet_edit.tif"

Here is my code:

## Plotting data points by coordinates (CH1903/LV03) on relief of switzerland

library(readxl)
library(raster)
library(ggplot2)

# import coordinates of weather stationes

coord <- read_excel("Q:/Projekte/Tarik/Coordinate_Plot/coordinate_plot/SMN-Stationen_20151222.xlsx")

# import raster (relief) file and modify it to be used in ggplot

rel <- raster("C:/Users/qualpra/Desktop/tif/landesgebiet_edit.tif")
rel_spdf <- as(rel, "SpatialPixelsDataFrame")

rel <- as.data.frame(rel_spdf)

head(rel)
landesgebiet_edit        x        y
               137 684552.0 295862.1
               122 684552.0 295777.5
               108 684636.7 295777.5
               116 684467.4 295692.8
               105 684552.0 295692.8
                93 684636.7 295692.8

ggplot() +
  geom_raster(data = rel, aes_string(x = "x", y = "y", alpha = "landesgebiet_edit")) +
  scale_alpha(name = "", range = c(0.6, 0), guide = F) +
  geom_point(data = coord, aes(x = X_KM_COORD_NU, y = Y_KM_COORD_NU), colour = "black", size = 3) +
  geom_point(data = coord, aes(x = X_KM_COORD_NU, y = Y_KM_COORD_NU), colour = "skyblue1", size = 1.5) +  
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) +
  theme(axis.title.y=element_blank(),
        axis.text.y=element_blank(),
        axis.ticks.y=element_blank()) +
  scale_fill_identity()

Output plot

Big thanks to Timo Grossenbacher!

score 0 · Answer 2 · answered Jan 07 '21 at 08:45

I came up with the same problem.

Here's my solution to edit the relief directly in R with rgdal:

library(ggplot2)
library(raster)
library(rgdal)

#Downloads from swisstopo:
r <-raster("..../20161001_SMV1000_SHAPE_CHLV95/Shapefiles/25_RELI.tif")

shapeCH <- readOGR("..../BOUNDARIES_2020_10/DATEN/swissBOUNDARIES3D/SHAPEFILE_LV03_LN02/swissBOUNDARIES3D_1_3_TLM_LANDESGEBIET.shp")
shapeCH <- spTransform(shapeCH, CRS(proj4string(r))) #Transform shape to raster CRS

#crop and mask
r2 <- crop(r, extent(shapeCH)) 
r3 <- rasterize(shapeCH, r2, mask=TRUE)

#ggplot relief
rel_spdf <- as(r3, "SpatialPixelsDataFrame")
rel <- as.data.frame(rel_spdf)

ggplot() +
  geom_raster(data = rel, aes_string(x = "x", y = "y", alpha = "X25_RELI"))

Add a relief (GEOtiff, .tif) to a ggplot of the swiss country borders ((polygon) shapefile, .shp)

2 Answers2