将浮雕 (GEOtiff, .tif) 添加到瑞士国家边界的 ggplot ((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 文件的形式拥有它。

当我尝试绘制光栅时出现错误："cannot allocate vector of size 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 在他的博客 post (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()

输出图

非常感谢 Timo Grossenbacher！

Answer 2

我遇到了同样的问题。

这是我使用 rgdal 在 R 中直接编辑浮雕的解决方案：

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"))

将浮雕 (GEOtiff, .tif) 添加到瑞士国家边界的 ggplot ((polygon) shapefile, .shp)

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

r

raster

shapefile

ggplot2

geotiff