定位Isoband/Isochrone内的点

huangapple
go评论108阅读模式
英文:

Locate points inside a Isoband/Isochrone

问题

以下是您要翻译的代码部分：

  1. map_500 <- ggplot()+
  2. geom_contour_filled(data = tapas_kriging_sf, aes(x = x, y = y, z=var1.pred), binwidth = 500)+ 
  3. geom_sf(data = world, fill = "transparent", color = "grey")+ 
  4. geom_sf(data = water, fill = "white")+ geom_sf(data = TAPAS_sf, size = 0.05)+ geom_sf(data = XYZ_sf, size = 0.05, color = "red")+ 
  5. coord_sf(xlim = c(XYZ_bbox$xmax+1,XYZ_bbox$xmin-1), ylim = c(XYZ_bbox$ymax+1, XYZ_bbox$ymin-1), expand = FALSE)+ 
  6. labs(title = "Isochrones Map of Europe and the Near East (TAPAS)", fill = "Isochrones") 
  7. map_500
  1. library(gstat)
  2. library(sf)
  3. library(readr)
  4. library(tidyverse)
  5. library(rnaturalearth)
  7. ##Data
  8. no2 <- read_csv(system.file("external/no2.csv", 
  9.                             package = "gstat"), show_col_types = FALSE)
  10. no2_sf <- st_as_sf(no2, crs = "OGC:CRS84", coords = c("station_longitude_deg", "station_latitude_deg")) 
  11. no2_sf <- st_transform(no2_sf, 32632)
  14. ##Kriging
  15. #Interpolation raster
  16. no2_bbox <- st_bbox(no2_sf)
  17. cell_size <- 10000
  18. x <- seq(no2_bbox$xmin, no2_bbox$xmax, by=cell_size)
  19. y <- seq(no2_bbox$ymin, no2_bbox$ymax, by=cell_size)
  20. no2_grid <- expand.grid(x=x, y=y)
  21. no2_grid$tmp <- 1
  22. plot(no2_grid$x, no2_grid$y, pch=19, cex=0.1)
  23. no2_grid <- st_as_sf(no2_grid, coords = c("x","y"), crs = st_crs(no2_sf))
  24. st_crs(no2_grid) <- st_crs(no2_sf)
  26. #Variogram
  27. no2_sample_variogram <- gstat::variogram(NO2~1, no2_sf)
  28. plot(no2_sample_variogram, plot.numbers = TRUE)
  29. no2_model_variogram <- vgm(psill = 16, "Exp", range = 200000, nugget = 1)
  30. plot(no2_sample_variogram, no2_model_variogram)
  31. no2_fit_variogram <- fit.variogram(no2_sample_variogram, no2_model_variogram)
  32. plot(no2_sample_variogram, no2_fit_variogram)
  34. #Ordinary Kriging
  35. no2_sf <- no2_sf[!duplicated(no2_sf$geometry),] #check which observation were removed
  36. no2_kriging <- gstat::krige(NO2~1, no2_sf, no2_grid, no2_fit_variogram)
  38. no2_kriging$x <- st_coordinates(no2_kriging)[,1]
  39. no2_kriging$y <- st_coordinates(no2_kriging)[,2]
  41. ##Ggplot
  42. 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))
  43. points <- st_as_sf(points, coords = c("x","y"), crs = st_crs(32632))
  45. germany  <- ne_countries(scale = "medium", returnclass = "sf", country = "Germany")
  46. germany <- st_transform(germany, 32632)
  48. ggplot()+
  49.   geom_contour_filled(data = no2_kriging, aes(x = x, y = y, z=var1.pred))+
  50.   geom_sf(data = germany,  fill = "transparent", color = "black")+
  51.   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 :

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

This produces the following image :

定位Isoband/Isochrone内的点

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.

  1.    var1.pred var1.var            geometry     x     y
  2. 1   6725.361       NA POINT (-10.5 34.79) -10.5 34.79
  3. 2   6730.202       NA POINT (-10.4 34.79) -10.4 34.79
  4. 3   6734.932       NA POINT (-10.3 34.79) -10.3 34.79
  5. 4   6739.531       NA POINT (-10.2 34.79) -10.2 34.79
  6. 5   6743.977       NA POINT (-10.1 34.79) -10.1 34.79
  7. 6   6748.251       NA POINT (-10   34.79) -10.0 34.79
  8. 7   6752.330       NA POINT (-9.9  34.79)  -9.9 34.79
  9. 8   6756.193       NA POINT (-9.8  34.79)  -9.8 34.79
  10. 9   6759.822       NA POINT (-9.7  34.79)  -9.7 34.79
  11. 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 定位Isoband/Isochrone内的点

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 :

  1. library(gstat)
  2. library(sf)
  3. library(readr)
  4. library(tidyverse)
  5. library(rnaturalearth)
  6. ##Data
  7. no2 &lt;- read_csv(system.file(&quot;external/no2.csv&quot;, 
  8. package = &quot;gstat&quot;), show_col_types = FALSE)
  9. no2_sf &lt;- st_as_sf(no2, crs = &quot;OGC:CRS84&quot;, coords = c(&quot;station_longitude_deg&quot;, &quot;station_latitude_deg&quot;)) 
  10. no2_sf &lt;- st_transform(no2_sf, 32632)
  11. ##Kriging
  12. #Interpolation raster
  13. no2_bbox &lt;- st_bbox(no2_sf)
  14. cell_size &lt;- 10000
  15. x &lt;- seq(no2_bbox$xmin, no2_bbox$xmax, by=cell_size)
  16. y &lt;- seq(no2_bbox$ymin, no2_bbox$ymax, by=cell_size)
  17. no2_grid &lt;- expand.grid(x=x, y=y)
  18. no2_grid$tmp &lt;- 1
  19. plot(no2_grid$x, no2_grid$y, pch=19, cex=0.1)
  20. no2_grid &lt;- st_as_sf(no2_grid, coords = c(&quot;x&quot;,&quot;y&quot;), crs = st_crs(no2_sf))
  21. st_crs(no2_grid) &lt;- st_crs(no2_sf)
  22. #Variogram
  23. no2_sample_variogram &lt;- gstat::variogram(NO2~1, no2_sf)
  24. plot(no2_sample_variogram, plot.numbers = TRUE)
  25. no2_model_variogram &lt;- vgm(psill = 16, &quot;Exp&quot;, range = 200000, nugget = 1)
  26. plot(no2_sample_variogram, no2_model_variogram)
  27. no2_fit_variogram &lt;- fit.variogram(no2_sample_variogram, no2_model_variogram)
  28. plot(no2_sample_variogram, no2_fit_variogram)
  29. #Ordinary Kriging
  30. no2_sf &lt;- no2_sf[!duplicated(no2_sf$geometry),] #check which observation were removed
  31. no2_kriging &lt;- gstat::krige(NO2~1, no2_sf, no2_grid, no2_fit_variogram)
  32. no2_kriging$x &lt;- st_coordinates(no2_kriging)[,1]
  33. no2_kriging$y &lt;- st_coordinates(no2_kriging)[,2]
  34. ##Ggplot
  35. points &lt;- data.frame(x = runif(10, min = no2_bbox$xmin, max = no2_bbox$xmax), y = runif(10, min = no2_bbox$ymin, max = no2_bbox$ymax))
  36. points &lt;- st_as_sf(points, coords = c(&quot;x&quot;,&quot;y&quot;), crs = st_crs(32632))
  37. germany  &lt;- ne_countries(scale = &quot;medium&quot;, returnclass = &quot;sf&quot;, country = &quot;Germany&quot;)
  38. germany &lt;- st_transform(germany, 32632)
  39. ggplot()+
  40. geom_contour_filled(data = no2_kriging, aes(x = x, y = y, z=var1.pred))+
  41. geom_sf(data = germany,  fill = &quot;transparent&quot;, color = &quot;black&quot;)+
  42. geom_sf(data = points, size = 0.5, color = &quot;red&quot;)

答案1

得分: 3

如果你简化一下流程，并使用stars栅格作为插值网格，然后在输出栅格上使用stars::st_contour()，你可以得到一个带有等时线多边形的sf对象。接着可以通过空间连接来对点进行分类。现在，在st_contour()调用中处理定义的箱子。

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

  1. library(sf)
  2. library(stars)
  3. library(gstat)
  4. library(rnaturalearth)
  5. library(readr)
  6. library(dplyr)
  7. library(ggplot2)
  9. ##数据
  10. no2 <- read_csv(system.file("external/no2.csv", 
  11.                             package = "gstat"), show_col_types = FALSE)
  12. no2_sf <- st_as_sf(no2, crs = "OGC:CRS84", coords = c("station_longitude_deg", "station_latitude_deg")) 
  13. no2_sf <- st_transform(no2_sf, 32632)
  15. ##克里金插值
  16. #插值栅格
  17. no2_grid <- st_bbox(no2_sf) %>% st_as_stars(dx = 10000)
  18. # stars对象:
  19. no2_grid 
  20. #> stars object with 2 dimensions and 1 attribute
  21. #> attribute(s):
  22. #>         Min. 1st Qu. Median Mean 3rd Qu. Max.
  23. #> values     0       0      0    0       0    0
  24. #> dimension(s):
  25. #>   from to  offset  delta                refsys x/y
  26. #> x    1 60  304638  10000 WGS 84 / UTM zone 32N [x]
  27. #> y    1 83 6086661 -10000 WGS 84 / UTM zone 32N [y]
  29. #半变异图
  30. no2_sample_variogram <- gstat::variogram(NO2~1, no2_sf)
  31. no2_model_variogram <- vgm(psill = 16, "Exp", range = 200000, nugget = 1)
  32. no2_fit_variogram <- fit.variogram(no2_sample_variogram, no2_model_variogram)
  34. #普通克里金
  35. # 现在返回 stars 栅格:
  36. no2_kriging <- gstat::krige(formula   = NO2~1, 
  37.                             locations = no2_sf, 
  38.                             newdata   = no2_grid, 
  39.                             model     = no2_fit_variogram)
  40. #> [using ordinary kriging]
  41. no2_kriging
  42. #> stars object with 2 dimensions and 2 attributes
  43. #> attribute(s):
  44. #>                 Min.  1st Qu.    Median      Mean   3rd Qu.     Max.
  45. #> var1.pred  2.6449424 7.281307  8.214847  8.569172  9.611575 19.09527
  46. #> var1.var   0.2019452 8.902331 11.587504 11.155864 13.871108 16.72409
  47. #> dimension(s):
  48. #>   from to  offset  delta                refsys x/y
  49. #> x    1 60  304638  10000 WGS 84 / UTM zone 32N [x]
  50. #> y    1 83 6086661 -10000 WGS 84 / UTM zone 32N [y]
  51. plot(no2_kriging)
  1. #从栅格中提取等时线，范围约为3到20，使用seq(0, 20, 5)来设置间隔
  2. isochrones <- st_contour(no2_kriging, breaks = seq(0, 20, 5))
  3. isochrones
  4. #> Simple feature collection with 4 features and 3 fields
  5. #> Geometry type: MULTIPOLYGON
  6. #> Dimension:     XY
  7. #> Bounding box:  xmin: 304638.2 ymin: 5256661 xmax: 904638.2 ymax: 6086661
  8. #> Projected CRS: WGS 84 / UTM zone 32N
  9. #>   var1.pred Min Max                       geometry
  10. #> 1     [0,5)   0   5 MULTIPOLYGON (((789638.2 60...
  11. #> 2    [5,10)   5  10 MULTIPOLYGON (((899638.2 52...
  12. #> 3   [10,15)  10  15 MULTIPOLYGON (((544131.7 58...
  13. #> 4   [15,20)  15  20 MULTIPOLYGON (((316428 5671...
  14. plot(isochrones["var1.pred"], key.pos = 1)
  1. ##采样点
  2. set.seed(123)
  4. points <- st_bbox(no2_grid) %>% st_as_sfc() %>% st_sample(100) %>% st_sf()
  6. # 空间连接以分类点，尝试从每个等时线中选择几个
  7. points <- st_join(points, isochrones["var1.pred"]) %>%
  8.   slice_head(n = 2, by = var1.pred)
  9. points
  10. #> Simple feature collection with 6 features and 1 field
  11. #> Geometry type: POINT
  12. #> Dimension:     XY
  13. #> Bounding box:  xmin: 319406.4 ymin: 5506681 xmax: 777621.3 ymax: 5754652
  14. #> Projected CRS: WGS 84 / UTM zone 32N
  15. #>   var1.pred                 geometry
  16. #> 1   [10,15) POINT (477184.7 5754652)
  17. #> 2   [10,15) POINT (576638.7 5506681)
  18. #> 3    [5,10) POINT (777621.3 5532905)
  19. #> 4    [5,10) POINT (550024.3 5662210)
  20. #> 5   [15,20) POINT (319406.4 5689204)
  21. #> 6   [15,20) POINT (371319.4 5739514)
  23. ## ggplot
  24. germany  <- ne_countries(scale = "medium", returnclass = "sf", country = "Germany")
  25. germany <- st_transform(germany, 32632)
  27. ggplot()+
  28.   geom_sf(data = isochrones, aes(fill =  var1.pred), alpha = .6) + 
  29.   geom_sf(data = germany,  fill = "transparent", color = "black")+
  30.   geom_sf(data = points, color = "red") +
  31.   geom_sf_label(data = points, aes(label = var1.pred), nudge_x = 60000, nudge_y = -25000, alpha =
  33. <details>
  34. <summary>英文:</summary>
  36. 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. 
  38. That extra column with matching isochrones for each point can be seen in the last code block.
  40. ``` r
  41. library(sf)
  42. library(stars)
  43. library(gstat)
  44. library(rnaturalearth)
  45. library(readr)
  46. library(dplyr)
  47. library(ggplot2)
  49. ##Data
  50. no2 &lt;- read_csv(system.file(&quot;external/no2.csv&quot;, 
  51.                             package = &quot;gstat&quot;), show_col_types = FALSE)
  52. no2_sf &lt;- st_as_sf(no2, crs = &quot;OGC:CRS84&quot;, coords = c(&quot;station_longitude_deg&quot;, &quot;station_latitude_deg&quot;)) 
  53. no2_sf &lt;- st_transform(no2_sf, 32632)
  55. ##Kriging
  56. #Interpolation raster
  57. no2_grid &lt;- st_bbox(no2_sf) %&gt;% st_as_stars(dx = 10000)
  58. # stars object:
  59. no2_grid 
  60. #&gt; stars object with 2 dimensions and 1 attribute
  61. #&gt; attribute(s):
  62. #&gt;         Min. 1st Qu. Median Mean 3rd Qu. Max.
  63. #&gt; values     0       0      0    0       0    0
  64. #&gt; dimension(s):
  65. #&gt;   from to  offset  delta                refsys x/y
  66. #&gt; x    1 60  304638  10000 WGS 84 / UTM zone 32N [x]
  67. #&gt; y    1 83 6086661 -10000 WGS 84 / UTM zone 32N [y]
  69. #Variogram
  70. no2_sample_variogram &lt;- gstat::variogram(NO2~1, no2_sf)
  71. no2_model_variogram &lt;- vgm(psill = 16, &quot;Exp&quot;, range = 200000, nugget = 1)
  72. no2_fit_variogram &lt;- fit.variogram(no2_sample_variogram, no2_model_variogram)
  74. #Ordinary Kriging 
  75. # now returns stars raster:
  76. no2_kriging &lt;- gstat::krige(formula   = NO2~1, 
  77.                             locations = no2_sf, 
  78.                             newdata   = no2_grid, 
  79.                             model     = no2_fit_variogram)
  80. #&gt; [using ordinary kriging]
  81. no2_kriging
  82. #&gt; stars object with 2 dimensions and 2 attributes
  83. #&gt; attribute(s):
  84. #&gt;                 Min.  1st Qu.    Median      Mean   3rd Qu.     Max.
  85. #&gt; var1.pred  2.6449424 7.281307  8.214847  8.569172  9.611575 19.09527
  86. #&gt; var1.var   0.2019452 8.902331 11.587504 11.155864 13.871108 16.72409
  87. #&gt; dimension(s):
  88. #&gt;   from to  offset  delta                refsys x/y
  89. #&gt; x    1 60  304638  10000 WGS 84 / UTM zone 32N [x]
  90. #&gt; y    1 83 6086661 -10000 WGS 84 / UTM zone 32N [y]
  91. plot(no2_kriging)

定位Isoband/Isochrone内的点<!-- -->

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

定位Isoband/Isochrone内的点<!-- -->

  1. ## Sample points
  2. set.seed(123)
  4. points &lt;- st_bbox(no2_grid) %&gt;% st_as_sfc() %&gt;% st_sample(100) %&gt;% st_sf()
  6. # spatial join to classify points, try to pick few from each isochrone
  7. points &lt;- st_join(points, isochrones[&quot;var1.pred&quot;]) %&gt;% 
  8.   slice_head(n = 2, by = var1.pred)
  9. points
  10. #&gt; Simple feature collection with 6 features and 1 field
  11. #&gt; Geometry type: POINT
  12. #&gt; Dimension:     XY
  13. #&gt; Bounding box:  xmin: 319406.4 ymin: 5506681 xmax: 777621.3 ymax: 5754652
  14. #&gt; Projected CRS: WGS 84 / UTM zone 32N
  15. #&gt;   var1.pred                 geometry
  16. #&gt; 1   [10,15) POINT (477184.7 5754652)
  17. #&gt; 2   [10,15) POINT (576638.7 5506681)
  18. #&gt; 3    [5,10) POINT (777621.3 5532905)
  19. #&gt; 4    [5,10) POINT (550024.3 5662210)
  20. #&gt; 5   [15,20) POINT (319406.4 5689204)
  21. #&gt; 6   [15,20) POINT (371319.4 5739514)
  23. ## ggplot
  24. germany  &lt;- ne_countries(scale = &quot;medium&quot;, returnclass = &quot;sf&quot;, country = &quot;Germany&quot;)
  25. germany &lt;- st_transform(germany, 32632)
  27. ggplot()+
  28.   geom_sf(data = isochrones, aes(fill =  var1.pred), alpha = .6) + 
  29.   geom_sf(data = germany,  fill = &quot;transparent&quot;, color = &quot;black&quot;)+
  30.   geom_sf(data = points, color = &quot;red&quot;) +
  31.   geom_sf_label(data = points, aes(label = var1.pred), nudge_x = 60000, nudge_y = -25000, alpha = .5) +
  32.   scale_fill_viridis_d() +
  33.   theme_minimal()

定位Isoband/Isochrone内的点<!-- -->

<sup>Created on 2023-07-11 with reprex v2.0.2</sup>

huangapple
  • 本文由 发表于 2023年7月11日 08:43:13
  • 转载请务必保留本文链接：https://go.coder-hub.com/76658081.html
  • geolocation
  • ggplot2
  • interpolation
  • r
  • r-sf
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定