在瑞士国家边界的ggplot中添加浮雕（GEOtiff，.tif）（（polygon）shapefile，.shp）

Question

我绘制了瑞士的区域（多边形shapefile）并通过坐标添加了点（瑞士气象站）。

# 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”下载数据。

在此地图上可以找到Countryborders（swissBOUNDARIES3D_1_3_TLM_KANTONSGEBIET.shp）：

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

在此地图上可以找到地形地形图（25_HYPSO.tif）：

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告诉我，他更改了CRS（从CH1903 + / LV95更改为CH1903 / LV03），并使用开源软件QGIS裁剪（导出了瑞士国家边界的范围）。

因此，这就是我对相同的.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()

输出图

非常感谢Timo Grossenbacher！