英文:
Locate points inside a Isoband/Isochrone
问题
以下是您要翻译的代码部分：
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
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")
希望这些翻译能够帮助您理解代码。
英文:
Part of my master's thesis involve doing some kriging spatial interpolation. I'm mapping the spread of agriculture into Europe from the Near east using an ordinary kriging interpolation over a map of Europe. I then rendered the map on ggplot. With the following code :
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
This produces the following image :
EDIT : Here's more information about the nature of the tapas_kriging_sf. It's an sf object, where each row is a point of the map with a value corresponding to the interpolation value (var1.pred) of that part of the map.
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
Now as you can see, the entirety of Europe is structured in colored "ribbons" called Isochrones or Isobands representing the arrival of agriculture at a give date (For example, in the yellow Isochrone, agriculture arrives at the earliest at the date of 8500 calBP...)
Scattered on this map, you can observe several red points. These points were not used to build the interpolation and come from another dataset.
My objective is : I'd like to create a column in this latter dataset which references each red points with the Isochrone it is located in.
For examples, for the red point in the Baleares (where you can see the arrow pointing), I'd like to have in the same row but in another column named "Isochrones" that it belongs to the Isochrone (7000,7500].
I figured out that my best bet would probably to polygonize the different isochrones and use st_intersect() or st_contains() to check for each point where they are located, but I have trouble creating a polygon of the Isochrones from the ggplot() map.
To better understand what I did, here's the entirety of the kriging process done with another data and semi-random point, as with ggplot.
Thanks in advance 
I used ggplot_build() to get the layer from the map containing the filled contour (Isochrone). I then try to convert it into an sf object and "polygonize" the different Isochrone but this led me nowhere.
EDIT :
I've reproduced a version of my code from start to finish with open data. The output is less perfect since my interpolation grid doesn't align specifically with my are of interest (Germany), in the original data I hide that by colouring the oceans in white to override the filled contour.
Here it is :
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")
答案1
得分: 3
如果你简化一下流程，并使用stars栅格作为插值网格，然后在输出栅格上使用stars::st_contour()，你可以得到一个带有等时线多边形的sf对象。接着可以通过空间连接来对点进行分类。现在，在st_contour()调用中处理定义的箱子。
对于每个点的匹配等时线的额外列可以在最后的代码块中看到。
library(sf)
library(stars)
library(gstat)
library(rnaturalearth)
library(readr)
library(dplyr)
library(ggplot2)
##数据
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)
##克里金插值
#插值栅格
no2_grid <- st_bbox(no2_sf) %>% st_as_stars(dx = 10000)
# stars对象:
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]
#半变异图
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)
#普通克里金
# 现在返回 stars 栅格:
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)
#从栅格中提取等时线，范围约为3到20，使用seq(0, 20, 5)来设置间隔
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)
##采样点
set.seed(123)
points <- st_bbox(no2_grid) %>% st_as_sfc() %>% st_sample(100) %>% st_sf()
# 空间连接以分类点，尝试从每个等时线中选择几个
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 =
<details>
<summary>英文:</summary>
If you simplify the process a bit and use `stars` raster as an interpolation grid, you can later use `stars::st_contour()` on the output raster to get `sf` object with isochrone polygons. Which in turn can be used to classify points through spatial join. Defining bins is now handled in `st_contour()` call.
That extra column with matching isochrones for each point can be seen in the last code block.
``` r
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()
<!-- -->
<sup>Created on 2023-07-11 with reprex v2.0.2</sup>
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论