在R中不规则边界的地图上绘制插值数据

Question

我在地图上绘制不规则边界的结果有困难。我使用IDW方法对气候数据进行了如下插值：

library(ggplot2)
library(gstat)
library(sp)
library(maptools)


df<-my_data

coordinates(df) = ~X + Y
x.range <- as.numeric(c(20.375, 31.375))
y.range <- as.numeric(c(52.375, 61.375))

grd2 <- expand.grid(X = seq(from = x.range[1], to = x.range[2], by = 0.01), Y = seq(from = y.range[1], to = y.range[2], by = 0.01))
coordinates(grd2) <- ~X + Y
gridded(grd2) <- TRUE
plot(grd2, cex = 1.5, col = "grey")
points(df, pch = 1, col = "red", cex = 1)

df2<-my_data
colnames(df2)<- c("ID", "lon", "lat", "vari")
idw2 <- idw(formula = vari ~ 1, locations = df, newdata = grd2)
idw.output2 = as.data.frame(idw2)
names(idw.output2)[1:3] <- c("lon", "lat", "var1.pred")

现在我想把这些idw的结果蒙上面具。在过去，我一直只使用一个国家的边界，没有任何问题。代码：

mymap <- readOGR("countryy.shp", layer="countryy")
summary(mymap)
wgsmap <- spTransform(mymap, CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
wgsmap.contour <- fortify(wgsmap)

library(raster)

idw.r <- rasterFromXYZ(idw.output2[, c("lon", "lat", "var1.pred")])
idw.crp <- crop(idw.r, wgsmap)
idw.msk <- mask(idw.crp, wgsmap)
idw.msk.dfr <- as.data.frame(rasterToPoints(idw.msk))
names(idw.msk.dfr)[1:2] <- c("lon", "lat")

#result
ggplot() + geom_tile(data = idw.msk.dfr, alpha = 0.8, aes(x = lon, y = lat, fill = round(var1.pred, 0))) + scale_fill_gradient(low = "cyan", high = "orange") + geom_path(data = wgsmap.contour, aes(long, lat, group = group), colour = "grey") + geom_point(data = df2, aes(x = lon, y = lat), shape = 18, colour = "red")

现在，我需要戴上一个新边界的面具，包括几个国家：

x.range <- as.numeric(c(20.375，31.375)) y.range <- as.numeric(c(52.375，61.375))

我不知道该怎么做。我是否应该根据我的坐标矩形来修剪整个世界的ESRI形状文件？或者我可以使用ggspatial的简单地图：

library("rnaturalearth")
library("rnaturalearthdata")
library(ggspatial)
world <- ne_countries(scale = "medium", returnclass = "sf")
class(world)

europe <- ggplot(data = world) +geom_sf() + coord_sf(xlim = c(20.375, 31.375), ylim = c(52.375, 61.375), expand = FALSE)

如果有人能帮我解决这个问题，我会非常感激的。

Answer 1

这个问题是用“GADMTools”软件包出售的。

joinmap <- gadm_sp_loadCountries(c("RUS","LTU","LVA", "EST", "FIN", "POL", "BLR"),  basefile = "./")

studyarea <- gadm_crop(joinmap, xmin=20.375, ymin=52.375, xmax=31.375, ymax=61.375)

gadm_exportToShapefile(studyarea, "path for shapefile")