Я собираюсь использовать library(mapdeck) для построения всего, в основном потому, что я разработал библиотеку, поэтому я очень хорошо с ней знаком. Он использует карты Mapbox, поэтому вам понадобится токен Mapbox.
Сначала получите данные
library(sf)
library(data.table)
fl <- sf::st_read("~/Documents/github/SO_geoSpatial_crop_Quest/GIS_cgbound/cgbound.shp") %>% sf::st_transform(crs = 4326)
gainsville_df <- fread("~/Documents/github/SO_geoSpatial_crop_Quest/new_dataframe_data.csv")
sf_gainsville <- sf::st_as_sf(gainsville_df, wkt = "Location")
## no need to transform, because it's already in Lon / Lat (?)
sf::st_crs( sf_gainsville ) <- 4326
#install.packages("tidycensus")
library(tidycensus)
tidycensus::census_api_key("21adc0b3d6e900378af9b7910d04110cdd38cd75", install = T, overwrite = T)
alachua <- tidycensus::get_acs(state = "FL", county = "Alachua", geography = "tract", geometry = T, variables = "B01003_001")
alachua <- sf::st_transform( alachua, crs = 4326 )
Это то, с чем мы работаем. Я строю полигоны и путь границы
library(mapdeck)
set_token( read.dcf("~/Documents/.googleAPI", fields = "MAPBOX"))
## this is what the polygons and the Alachua boundary looks like
mapdeck() %>%
add_polygon(
data = alachua
, fill_colour = "NAME"
) %>%
add_path(
data = fl
, stroke_width = 50
)
Для начала я собираюсь сделать полигон границы
boundary_poly <- sf::st_cast(fl, "POLYGON")
Тогда мы можем получить эти многоугольники полностью внутри границы
idx <- sf::st_contains(
x = boundary_poly
, y = alachua
)
idx <- unlist( sapply( idx, `[`) )
sf_contain <- alachua[ idx, ]
mapdeck() %>%
add_polygon(
data = sf_contain
, fill_colour = "NAME"
) %>%
add_path(
data = fl
)
И те, которые «касаются» границы
idx <- sf::st_crosses(
x = fl
, y = alachua
)
idx <- unlist( idx )
sf_crosses <- alachua[ idx, ]
mapdeck() %>%
add_polygon(
data = sf_crosses
, fill_colour = "NAME"
) %>%
add_path(
data = fl
)
Те, которые полностью снаружи - это многоугольники, которые не касаются границы и не находятся внутри нее
sf_outside <- sf::st_difference(
x = alachua
, y = sf::st_union( sf_crosses )
)
sf_outside <- sf::st_difference(
x = sf_outside
, y= sf::st_union( sf_contain )
)
mapdeck() %>%
add_polygon(
data = sf_outside
, fill_colour = "NAME"
) %>%
add_path(
data = fl
)
нам нужен способ «вырезать» те, которые касаются границы (sf_crosses), чтобы у нас был «внутренний» и «внешний» раздел для каждого многоугольника
Нам нужно работать с каждым многоугольником за раз и «разбивать» его по линиям, которые пересекают его.
Может быть способ сделать это с помощью lwgeom::st_split, но я продолжал получать ошибки
Чтобы помочь с этим, я использую версию для разработки моей sfheaders библиотеки
# devtools::install_github("dcooley/sfheaders")
res <- lapply( 1:nrow( sf_crosses ), function(x) {
## get the intersection of the polygon and the boundary
sf_int <- sf::st_intersection(
x = sf_crosses[x, ]
, y = fl
)
## we only need lines, not MULTILINES
sf_lines <- sfheaders::sf_cast(
sf_int, "LINESTRING"
)
## put a small buffer around the lines to make them polygons
sf_polys <- sf::st_buffer( sf_lines, dist = 0.0005 )
## Find the difference of these buffers and the polygon
sf_diff <- sf::st_difference(
sf_crosses[x, ]
, sf::st_union( sf_polys )
)
## this result is a MULTIPOLYGON, which is the original polygon from
## sf_crosses[x, ], split by the lines which cross it
sf_diff
})
## The result of this is all the polygons which touch the boundary path have been split
sf_res <- do.call(rbind, res)
, поэтому sf_res теперь должны быть всеми полигонами которые «касаются» пути, но разделяются там, где путь пересекает их
mapdeck() %>%
add_polygon(
data = sf_res
, stroke_colour = "#FFFFFF"
, stroke_width = 100
) %>%
add_path(
data = fl
, stroke_colour = "#FF00FF"
)
И мы можем увидеть это, увеличив
Теперь мы можем найти, какие из них находятся внутри и снаружи пути
sf_in <- sf::st_join(
x = sf_res
, y = boundary_poly
, left = FALSE
)
sf_out <- sf::st_difference(
x = sf_res
, y = sf::st_union( boundary_poly )
)
mapdeck() %>%
add_path(
data = fl
, stroke_width = 50
, stroke_colour = "#000000"
) %>%
add_polygon(
data = sf_in
, fill_colour = "NAME"
, palette = "viridis"
, layer_id = "in"
) %>%
add_polygon(
data = sf_out
, fill_colour = "NAME"
, palette = "plasma"
, layer_id = "out"
)
Теперь у нас есть все объекты, которые нас интересуют
sf_contain - все многоугольники полностью внутри бондара sf_in - все многоугольники касаются границы внутри sf_out - все многоугольники, касающиеся границы снаружи sf_outside - все остальные многоугольники
mapdeck() %>%
add_path(
data = fl
, stroke_width = 50
, stroke_colour = "#000000"
) %>%
add_polygon(
data = sf_contain
, fill_colour = "NAME"
, palette = "viridis"
, layer_id = "contained_within_boundary"
) %>%
add_polygon(
data = sf_in
, fill_colour = "NAME"
, palette = "cividis"
, layer_id = "touching_boundary_inside"
) %>%
add_polygon(
data = sf_out
, fill_colour = "NAME"
, palette = "plasma"
, layer_id = "touching_boundary_outside"
) %>%
add_polygon(
data = sf_outside
, fill_colour = "NAME"
, palette = "viridis"
, layer_id = "outside_boundary"
)
