Блог

Здравствуйте (здесь новичок в картографировании!)

У меня была хорошая охота, и я не могу найти решение того, что кажется сложной проблемой.

У меня есть базовые данные XY (координаты):

Что я хочу сделать, так это создать соседние полигоны на основе данных координат, которые не перекрываются и имеют определенный предел размера (чтобы они не простирались навсегда в океан).

Простите мои плохие навыки MS Paint, но желаемый результат был бы примерно таким:

У меня есть полигон, который отмечает границу раздела земля / море, поэтому полигоны также не могут перекрывать его.

Я использую листовку, чтобы сделать эти карты интерактивными, это не для какого-либо статистического анализа, а для предоставления обзора.

Конечная цель состоит в том, чтобы каждый полигон был окрашен переменной (например, температурой) с наложением экологических данных.

Некоторые примеры данных:

     > 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')
  

Меня весь день расстраивало, что в большинстве примеров используются загруженные полигоны (например, то, что кажется классическими штатами США, а не их создание)

Любая помощь приветствуется! (или я прошу слишком многого!)

Джим

Вот кое-что, с чего можно начать:

 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 ) 
  

Приведенный выше ответ мне немного помог, но я нашел свое решение здесь: 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