r - 将浮雕(GEOtiff，.tif)添加到瑞士国家边界的 ggplot((多边形)shapefile，.shp)

Question

我绘制了瑞士的面积(多边形形状文件)并通过坐标添加了点(瑞士气象站)。

# 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()

现在我想在图形下面添加地形浮雕，我拥有 GEOtiff 文件形式的地形浮雕。

当我尝试绘制栅格时，出现错误:“无法分配大小为 196.7 Mb 的向量”。

另一个问题是，这两个文件的范围不同。所以我不能按程度裁剪它们。有没有办法通过他们的 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"  

您可以从联邦地形局“swisstopo”下载数据。

国家边界 (swissBOUNDARIES3D_1_3_TLM_KANTONSGEBIET.shp) 可以在此 map 中找到:

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

地形起伏 (25_HYPSO.tif) 可以在此 map 中找到:

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

最终产品应该看起来像 Timo Grossenbacher 在他的博客文章 ( https://timogrossenbacher.ch/2016/12/beautiful-thematic-maps-with-ggplot2-only/ ) 中所做的那样。我尝试了他的方法，但没有成功，因为他提前编辑了他的救济文件，而我不知道如何编辑。

Answer 1

我现在成功了。 Timo Grossenbacher 告诉我，他使用开源软件 QGIS 更改了 CRS(从 CH1903+/LV95 到 CH1903/LV03)并进行了裁剪(导出瑞士国家边界的范围)。

这就是我对同一个 .tif 文件和另一个用于导出的 .shp 文件所做的操作，该文件仅包含国家/地区边界(/swissBOUNDARIES3D_1_3_TLM_LANDESGEBIET.shp，可以在同一下载中找到)。

我将生成的文件命名为:“landesgebiet_edit.tif”

这是我的代码:

## 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()

输出图

非常感谢蒂莫·格罗森巴赫!