从空间点数据创建边界多边形以进行绘图

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 Earth 上手动找到

## 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.

创建 Voroni 图

## 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