从空间点数据创建边界多边形以在传单中绘制

Question

你好（这里是地图新手！）

我已经进行了很好的搜索，但找不到解决这个看似棘手的问题的方法。

我有基本的 XY（坐标）数据：

我想要做的是根据不重叠且具有一定大小限制的坐标数据创建相邻的多边形（因此它们不会永远延伸到海洋中）。

请原谅我糟糕的 MS Paint 技能，但理想的结果应该是这样的：

我有一个标记陆地/海洋界面的多边形，因此多边形也不能重叠。

我正在使用 Leaflet 使这些地图具有交互性，它不是用于任何统计分析，而是提供概述。

最终目标是让每个多边形由一个变量（例如温度）着色，并覆盖生态数据。

一些示例数据：

    > data[1:10,]
   Station  Lat_dec  Long_dec Surface_T
1      247 50.33445 -2.240283     15.19
2      245 50.58483 -2.535217     14.11
3      239 50.16883 -2.509250     15.41
4      225 50.32848 -2.765967     15.34
5      229 50.63900 -2.964800     14.09
6      227 50.33757 -3.303217     15.12
7      217 50.16657 -3.563817     15.13
8      207 49.66683 -3.556550     15.04
9      213 50.16512 -3.824667     14.97
10     219 49.83707 -3.815483     14.78

生成图 1 的代码是一个基本的传单脚本：

leaflet() %>% 
  addProviderTiles('Esri.OceanBasemap'
  ) %>% 
  addCircleMarkers(data = data,
                   lng = ~Long_dec,
                   lat = ~Lat_dec,
                   radius = 2
  ) %>%
  addPolygons(data = Land,
              weight = 1,
              color = 'black')

整天让我很沮丧，大多数示例都使用下载的多边形（例如，似乎是经典的美国州，而不是制作它们）

非常感谢任何帮助！ （还是我要求太多了！）

吉姆

Answer 1

这里有一些东西可以开始：

library(sf)
library(dplyr)

#create sf object with points
stations <- st_as_sf( df, coords = c( "Long_dec", "Lat_dec" ) ) 

#create voronoi/thiessen polygons
v <- stations %>% 
  st_union() %>%
  st_voronoi() %>%
  st_collection_extract()

library(leaflet)
leaflet() %>% 
  addTiles() %>% 
  addCircleMarkers( data = stations ) %>%
  addPolygons( data = v ) 

Answer 2

上述答案对我有所帮助，但我在这里找到了我的解决方案：https://github.com/r-spatial/sf/issues/474，并认为我会分享我正在研究的示例，以防它帮助其他人并使用可能感兴趣的真实数据。

## reshape used for colsplit function only
library(reshape2)
library(dplyr)
library(sf)

##define coordinate systems (rdnew is for Netherlands, but works for Denmark also)
latlong = "+init=epsg:4326"
rdnew = "+init=epsg:28992"

这些向量创建了来自 solarheatdata.eu 的太阳能热区供热站点的数据框，但最终必须在 Google 地球上手动找到

## create vector of names
site_names <- c("Strandby","Taars","Vraa","Jerslev", "Saeby", "Saeby 2","Jetsmark","Aabybro","Hjallerup", 
                "Dronninglund","Asaa","Ulsted", "Mou", "Snedsted", "Durup","Hadsund","Gjerlev","Ramsing-Lem-Lihme", "Haderup",
                "Feldborg","Frederiks","Karup" ,"Silkeborg", "Rye","Braedstrup 2","Braedstrup","Ejstrupholm", 
                "Torring" , "Torring2","Vildbjerg","Ornhoj-Gronbjerg", "Tim", "Ringkjobing", "Egtved" ,"Hejnsvig", 
                "Skovlund", "Tistrup" ,"Sig" , "Oksbol" , "Gording", "Holsted" , "Christiansfeld", "Gram" ,"Vojens",
                "Toftlund","Logumkloster", "Grasten" ,"Broager", "St. Rise","Aeroskobing","Marstel","Sydlangeland",
                "Tommerup","Ringe", "Lolland","Oster", "Sydfalster","Refa-Gedser","Stege", "Sandved-Tornemark", "Fuglebjerg", 
                "Haslev" ,"Hvidebaek", "Svebolle","Smorum", "Skuldelev","Jaegerspris","Nykobings", "Hundested", "Vejby", 
                "Helsinge")

## Vector of IDs

site_ids <- c("4", "79","45","51","14","115","49","110","81","37","40","3","32","47","114","54","53","112","44","24","30", 
"31", "100", "43", "99", "1", "16",  "5",  "86", "42", "25", "27", "11", "76", "12", "23", "15", "34", "13", "19",  "98",  "28", 
"6", "18", "36","50", "22", "10", "39", "8", "2",  "33", "82","116", "97", "93","17", "102", "83", "94", "96","103", "35",
"29", "104", "48", "9",  "41",  "46","21", "20") 

## site coordinates
site_latlons <- c("57.488519, 10.486884", "57.392246, 10.120313","57.361012, 9.947647",  "57.284072, 10.096013", "57.317312, 10.502039",
                  "57.315924, 10.501388", "57.204279, 9.670758",  "57.146419, 9.762861",  "57.175564, 10.151121", "57.162937, 10.253277",
                  "57.155180, 10.394182", "57.078654, 10.252327", "56.962971, 10.214036", "56.888538, 8.538094",  "56.741672, 8.963102",
                  "56.736111, 10.112385", "56.591002, 10.136062", "56.591399, 8.791363", "56.392732, 8.990445",  "56.336781, 8.940749",
                  "56.337432, 9.245461",  "56.311003, 9.163035",  "56.208234, 9.546451",  "56.069370, 9.692441", "55.977873, 9.623761",
                  "55.978695, 9.629253",  "55.980475, 9.293377",  "55.855816, 9.492632",  "55.854990, 9.496065",  "56.206777, 8.773864",
                  "56.204275, 8.573454",  "56.190090, 8.312271",  "56.094087, 8.284403",  "55.620152, 9.323732",  "55.695530, 9.007244",
                  "55.742032, 8.702813", "55.717618, 8.613441",  "55.667168, 8.565599",  "55.619889, 8.288348",  "55.476738, 8.803819",
                  "55.506467, 8.903544",  "55.368032, 9.488627", "55.276828, 9.049574",  "55.241101, 9.285634",  "55.199975, 9.080764",
                  "55.051287, 8.960624",  "54.932625, 9.604330",  "54.884324, 9.676783", "54.852033, 10.408747", "54.883866, 10.415010",
                  "54.852505, 10.506659", "54.798816, 10.710182", "55.329109, 10.199913", "55.252603, 10.496850", "54.825794, 11.171034",
                  "54.758955, 11.818087", "54.718146, 11.923282", "54.577268, 11.935755", "55.000648, 12.291542", "55.254072, 11.522817",
                  "55.298998, 11.535034", "55.333998, 11.970144", "55.613226, 11.209546", "55.654732, 11.292038", "55.729510, 12.307020",
                  "55.775250, 12.016991", "55.829130, 12.000715", "55.917183, 11.651009", "55.963456, 11.880701", "56.069852, 12.152346",
                  "56.009482, 12.205189")

##combine and create spatial data frame
sites_sf <- data.frame(site = site_names, id = site_ids, latlon = site_latlons) %>% 
  mutate(colsplit(latlon, ",", c("lon", "lat"))) %>%
  st_as_sf(coords = c("lat", "lon"), crs = 4326) %>% 
  st_transform(28992) %>% 
  dplyr::select(-latlon)

rnaturalearth 和 rnatural earthdata 包可以轻松导入世界各国的形状文件。高分辨率可以下载，但需要安装 rnaturalearthhires，第一次经常失败，所以这里选择了“中”。

## shapefile of Denmark
library(rnaturalearth)

## download medium scale sf polygon of denmark. 
dk <- ne_countries(country = 'denmark', scale = 'medium', returnclass = 'sf')

## convert to a metres crs
dk_rd <- dk %>% 
  st_transform(28992) %>% 
  st_buffer(5000) ## as the shapefile is not that high resolution some edges are missed off. Buffer is used to cover more coastline, but is not essential.

创建沃罗尼图

## combine sites
d <- st_union(sites_sf)
## create voroni
v <- st_voronoi(d)

## trim the plot to the country shapefile
p1 <- st_as_sf(st_intersection(st_cast(v), dk_rd), crs = 28992)

## tell R that the geometry is polygons
p1 <- st_cast(p1, "MULTIPOLYGON")

plot(p1)
plot(sites_sf, add = TRUE)



## see on a map
library(mapview)

mapview(p1)+sites_sf