Get relief of Switzerland from a shp-file-CodePudding

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

CodePudding user response：

A simple way to get 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)

You can do the same things with the Swiss government data from the website you point to in your 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)

CodePudding user response：

Thanks to your comments - it works!

Downloading the asc-dataset from https://www.swisstopo.admin.ch/en/geodata/height/dhm25200.html and

relief <- raster(DHM200.asc) %>%  
  as("SpatialPixelsDataFrame") %>%
  as.data.frame()

ggplot()  
  
  # first: draw the relief
  geom_raster(
    data = relief,
    inherit.aes = FALSE,
    aes(
      x = x,
      y = y,
      alpha = DHM200
    )
  )

gives the desired map!