r - 定位等值带/等时线内的点

Question

我的硕士论文的一部分涉及进行一些克里格空间插值。我正在使用欧洲 map 上的普通克里金插值法来绘制农业从近东到欧洲的传播情况。然后我在 ggplot 上渲染了 map 。使用以下代码:

map_500 <- ggplot()+
geom_contour_filled(data = tapas_kriging_sf, aes(x = x, y = y, z=var1.pred), binwidth = 500)+ 
geom_sf(data = world, fill = "transparent", color = "grey")+ 
geom_sf(data = water, fill = "white")+ geom_sf(data = TAPAS_sf, size = 0.05)+ geom_sf(data = XYZ_sf, size = 0.05, color = "red")+ 
coord_sf(xlim = c(XYZ_bbox$xmax+1,XYZ_bbox$xmin-1), ylim = c(XYZ_bbox$ymax+1, XYZ_bbox$ymin-1), expand = FALSE)+ 
labs(title = "Isochrones Map of Europe and the Near East (TAPAS)", fill = "Isochrones") 
map_500

这会产生以下图像:

编辑:以下是有关 Tapas_kriging_sf 性质的更多信息。它是一个 sf 对象，其中每一行都是 map 上的一个点，其值对应于 map 该部分的插值值 (var1.pred)。

   var1.pred var1.var            geometry     x     y
1   6725.361       NA POINT (-10.5 34.79) -10.5 34.79
2   6730.202       NA POINT (-10.4 34.79) -10.4 34.79
3   6734.932       NA POINT (-10.3 34.79) -10.3 34.79
4   6739.531       NA POINT (-10.2 34.79) -10.2 34.79
5   6743.977       NA POINT (-10.1 34.79) -10.1 34.79
6   6748.251       NA POINT (-10   34.79) -10.0 34.79
7   6752.330       NA POINT (-9.9  34.79)  -9.9 34.79
8   6756.193       NA POINT (-9.8  34.79)  -9.8 34.79
9   6759.822       NA POINT (-9.7  34.79)  -9.7 34.79
10  6763.195       NA POINT (-9.6  34.79)  -9.6 34.79

现在，如您所见，整个欧洲由称为等时线或等带线的彩色“丝带”构成，代表农业在给定日期的到达(例如，在黄色等时线中，农业最早到达的日期8500 calBP...)

在这张 map 上，你可以观察到几个红​​点。这些点未用于构建插值，而是来自另一个数据集。

我的目标是:我想在后一个数据集中创建一个列，该列引用每个红点及其所在的等时线。

例如，对于巴利阿里群岛中的红点(您可以看到箭头指向的位置)，我希望在同一行但在另一列名为“等时线”的列中，它属于等时线(7000， 7500]。

我发现我最好的选择可能是多边形化不同的等时线并使用 st_intersect() 或 st_contains() 来检查它们所在的每个点，但我无法从 ggplot( 创建等时线的多边形) ) map 。

为了更好地理解我所做的，这里是使用另一个数据和半随机点完成的整个克里金过程，就像 ggplot 一样。

提前致谢:)

我使用 ggplot_build() 从包含填充轮廓(等时线)的 map 中获取图层。然后我尝试将其转换为 sf 对象并“多边形化”不同的等时线，但这使我无处可去。

编辑: 我用开放数据从头到尾重现了我的代码版本。输出不太完美，因为我的插值网格与我感兴趣的区域(德国)没有特别对齐，在原始数据中，我通过将海洋着色为白色以覆盖填充的轮廓来隐藏它。 这是:

library(gstat)
library(sf)
library(readr)
library(tidyverse)
library(rnaturalearth)

##Data
no2 <- read_csv(system.file("external/no2.csv", 
                            package = "gstat"), show_col_types = FALSE)
no2_sf <- st_as_sf(no2, crs = "OGC:CRS84", coords = c("station_longitude_deg", "station_latitude_deg")) 
no2_sf <- st_transform(no2_sf, 32632)


##Kriging
#Interpolation raster
no2_bbox <- st_bbox(no2_sf)
cell_size <- 10000
x <- seq(no2_bbox$xmin, no2_bbox$xmax, by=cell_size)
y <- seq(no2_bbox$ymin, no2_bbox$ymax, by=cell_size)
no2_grid <- expand.grid(x=x, y=y)
no2_grid$tmp <- 1
plot(no2_grid$x, no2_grid$y, pch=19, cex=0.1)
no2_grid <- st_as_sf(no2_grid, coords = c("x","y"), crs = st_crs(no2_sf))
st_crs(no2_grid) <- st_crs(no2_sf)

#Variogram
no2_sample_variogram <- gstat::variogram(NO2~1, no2_sf)
plot(no2_sample_variogram, plot.numbers = TRUE)
no2_model_variogram <- vgm(psill = 16, "Exp", range = 200000, nugget = 1)
plot(no2_sample_variogram, no2_model_variogram)
no2_fit_variogram <- fit.variogram(no2_sample_variogram, no2_model_variogram)
plot(no2_sample_variogram, no2_fit_variogram)

#Ordinary Kriging
no2_sf <- no2_sf[!duplicated(no2_sf$geometry),] #check which observation were removed
no2_kriging <- gstat::krige(NO2~1, no2_sf, no2_grid, no2_fit_variogram)

no2_kriging$x <- st_coordinates(no2_kriging)[,1]
no2_kriging$y <- st_coordinates(no2_kriging)[,2]

##Ggplot
points <- data.frame(x = runif(10, min = no2_bbox$xmin, max = no2_bbox$xmax), y = runif(10, min = no2_bbox$ymin, max = no2_bbox$ymax))
points <- st_as_sf(points, coords = c("x","y"), crs = st_crs(32632))
  
germany  <- ne_countries(scale = "medium", returnclass = "sf", country = "Germany")
germany <- st_transform(germany, 32632)

ggplot()+
  geom_contour_filled(data = no2_kriging, aes(x = x, y = y, z=var1.pred))+
  geom_sf(data = germany,  fill = "transparent", color = "black")+
  geom_sf(data = points, size = 0.5, color = "red")

Answer 1

如果您稍微简化一下过程并使用 stars栅格作为插值网格，您可以稍后使用 stars::st_contour()在输出栅格上获取 sf具有等时线多边形的对象。反过来又可用于通过空间连接对点进行分类。定义垃圾箱现在在 st_contour() 中处理打电话。

在最后一个代码块中可以看到每个点具有匹配等时线的额外列。

library(sf)
library(stars)
library(gstat)
library(rnaturalearth)
library(readr)
library(dplyr)
library(ggplot2)

##Data
no2 <- read_csv(system.file("external/no2.csv", 
                            package = "gstat"), show_col_types = FALSE)
no2_sf <- st_as_sf(no2, crs = "OGC:CRS84", coords = c("station_longitude_deg", "station_latitude_deg")) 
no2_sf <- st_transform(no2_sf, 32632)

##Kriging
#Interpolation raster
no2_grid <- st_bbox(no2_sf) %>% st_as_stars(dx = 10000)
# stars object:
no2_grid 
#> stars object with 2 dimensions and 1 attribute
#> attribute(s):
#>         Min. 1st Qu. Median Mean 3rd Qu. Max.
#> values     0       0      0    0       0    0
#> dimension(s):
#>   from to  offset  delta                refsys x/y
#> x    1 60  304638  10000 WGS 84 / UTM zone 32N [x]
#> y    1 83 6086661 -10000 WGS 84 / UTM zone 32N [y]

#Variogram
no2_sample_variogram <- gstat::variogram(NO2~1, no2_sf)
no2_model_variogram <- vgm(psill = 16, "Exp", range = 200000, nugget = 1)
no2_fit_variogram <- fit.variogram(no2_sample_variogram, no2_model_variogram)

#Ordinary Kriging 
# now returns stars raster:
no2_kriging <- gstat::krige(formula   = NO2~1, 
                            locations = no2_sf, 
                            newdata   = no2_grid, 
                            model     = no2_fit_variogram)
#> [using ordinary kriging]
no2_kriging
#> stars object with 2 dimensions and 2 attributes
#> attribute(s):
#>                 Min.  1st Qu.    Median      Mean   3rd Qu.     Max.
#> var1.pred  2.6449424 7.281307  8.214847  8.569172  9.611575 19.09527
#> var1.var   0.2019452 8.902331 11.587504 11.155864 13.871108 16.72409
#> dimension(s):
#>   from to  offset  delta                refsys x/y
#> x    1 60  304638  10000 WGS 84 / UTM zone 32N [x]
#> y    1 83 6086661 -10000 WGS 84 / UTM zone 32N [y]
plot(no2_kriging)

# contrours from raster, range is ~ 3 .. 20, using seq(0, 20, 5) for breaks
isochrones <- st_contour(no2_kriging, breaks = seq(0, 20, 5))
isochrones
#> Simple feature collection with 4 features and 3 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 304638.2 ymin: 5256661 xmax: 904638.2 ymax: 6086661
#> Projected CRS: WGS 84 / UTM zone 32N
#>   var1.pred Min Max                       geometry
#> 1     [0,5)   0   5 MULTIPOLYGON (((789638.2 60...
#> 2    [5,10)   5  10 MULTIPOLYGON (((899638.2 52...
#> 3   [10,15)  10  15 MULTIPOLYGON (((544131.7 58...
#> 4   [15,20)  15  20 MULTIPOLYGON (((316428 5671...
plot(isochrones["var1.pred"], key.pos = 1)

## Sample points
set.seed(123)

points <- st_bbox(no2_grid) %>% st_as_sfc() %>% st_sample(100) %>% st_sf()

# spatial join to classify points, try to pick few from each isochrone
points <- st_join(points, isochrones["var1.pred"]) %>% 
  slice_head(n = 2, by = var1.pred)
points
#> Simple feature collection with 6 features and 1 field
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 319406.4 ymin: 5506681 xmax: 777621.3 ymax: 5754652
#> Projected CRS: WGS 84 / UTM zone 32N
#>   var1.pred                 geometry
#> 1   [10,15) POINT (477184.7 5754652)
#> 2   [10,15) POINT (576638.7 5506681)
#> 3    [5,10) POINT (777621.3 5532905)
#> 4    [5,10) POINT (550024.3 5662210)
#> 5   [15,20) POINT (319406.4 5689204)
#> 6   [15,20) POINT (371319.4 5739514)

## ggplot
germany  <- ne_countries(scale = "medium", returnclass = "sf", country = "Germany")
germany <- st_transform(germany, 32632)

ggplot()+
  geom_sf(data = isochrones, aes(fill =  var1.pred), alpha = .6) + 
  geom_sf(data = germany,  fill = "transparent", color = "black")+
  geom_sf(data = points, color = "red") +
  geom_sf_label(data = points, aes(label = var1.pred), nudge_x = 60000, nudge_y = -25000, alpha = .5) +
  scale_fill_viridis_d() +
  theme_minimal()

^{创建于 2023 年 7 月 11 日 reprex v2.0.2}