我绘制了瑞士的区域(多边形形状文件)并通过坐标添加了点(瑞士气象站)。
# 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) 可以在这张 map 中找到:
https://shop.swisstopo.admin.ch/de/products/landscape/boundaries3D
可以在此 map 中找到地形起伏 (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/ ) 中所做的那样。我尝试了他的方法，但没有成功，因为他提前编辑了他的救济文件，我不知道怎么做。
最佳答案
我现在设法做到了。 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!
关于r - 将浮雕(GEOtiff，.tif)添加到瑞士国家边界的ggplot((多边形)shapefile，.shp)，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/49557445/