Search code examples
rgisrastershapefile

Get relief of Switzerland from a shp-file

I would like an official source for drawing the relief over a map of Switzerland.

I downloaded https://cms.geo.admin.ch/ogd/topography/DHM25_BM_SHP.zip to use the dhm25_p.shp-File

Now using the code

aux <- st_read(dsn="dhm25_p.shp")
auxx <- as_Spatial(aux)
auxxx <- as.data.frame(auxx)

ggplot() +
  
  # draw the relief
  geom_raster(
    data = auxxx,
    inherit.aes = FALSE,
    aes(
      x = coords.x1,
      y = coords.x2,
      alpha = coords.x3
    )
  )

I'll get the error

Error in `geom_raster()`:
! Problem while converting geom to grob.
ℹ Error occurred in the 1st layer.
Caused by error:
! cannot allocate vector of size 6539542.5 Gb
Solution

  • A shapefile of points is generally not a good source for mapping continuous phenomena such as elevation. It is better to use raster data.

    Here is a simple way to make an elevation map for Switzerland

    library(terra)
library(geodata)
x <- geodata::elevation_30s("Switzerland", ".")
plot(x)

    Add contour lines

    v <- as.contour(x, levels=c(500,1000,3000))
lines(v)

    Or show shaded relief

    slope <- terrain(x, "slope", unit="radians")
aspect <- terrain(x, "aspect", unit="radians")
hill <- shade(slope, aspect, 40, 270)
plot(hill, col=gray(seq(0,1,.01)), legend=F, axes=F, mar=1)

    To show relief and elevation

    plot(hill, col=gray(seq(0,1,.01)), legend=F, axes=F)
plot(x, col=terrain.colors(25, alpha=.5), add=T, axes=F, legend=T)

    enter image description here

    You can do the same things with the Swiss government data from the website you point to in your (now hidden) answer. But it takes some more work if you want to remove the areas outside of Switzerland.

    y <- rast("DHM200.asc")
# Assign the coordinate reference system (Landesvermessung 1903)
crs(y) <- "EPSG:21781" 

# get the outline of the country and project it to the crs of the raster data 
swz <- geodata::gadm("Switzerland", level=1, path=".")
pswz <- project(swz, y)

# remove values for areas outside of Switzerland
y <- mask(y, pswz)

plot(y)
lines (pswz)

    And with ggplot you can use geom_spatraster

    library(tidyterra)
library(ggplot2)
ggplot() + geom_spatraster(data = y)