如何保留空间特征的唯一交集并删除边界外的所有内容？

Question

我试图摆脱落在我读取的 shapefile 边界之外的空间几何。如果没有像 Photoshop 这样的手动软件，是否可以做到这一点？或者我手动移除跨越城市边界之外的区域。例如，我拿出了 14 张小册子，这是有结果的：

我已经提供了数据的所有子集和自己测试的密钥。代码脚本如下，数据集为https://github.com/THsTestingGround/SO_geoSpatial_crop_Quest。

我st_intersection(gainsville_df$Geomtry$x, gnv_poly$geometry)在转换Geomtry为sf.

library(sf)
library(tigris)
library(tidyverse)
library(tidycensus)
library(readr)
library(data.table)

#reading the shapefile
gnv_poly <-  sf::st_read("PATH\\GIS_cgbound\\cgbound.shp") %>% 
                sf::st_transform(crs = 4326) %>% 
                sf::st_polygonize() %>% 
                sf::st_union()

#I have taken the "geometry" of latitude and longitude because it was corrupting my csv, but we can rebuild like so
gnv_latlon <- readr::read_csv("new_dataframe_data.csv") %>% 
                dplyr::select(ID,
                              Latitude,
                              Longitude,
                              Location) %>%
                dplyr::mutate(Location = gsub(x= Location, pattern = "POINT \\(|\\)", replacement = "")) %>% 
                tidyr::separate(col = "Location", into = c("lon", "lat"), sep = " ") %>% 
                sf::st_as_sf(coords = c(4,5)) %>% 
                sf::st_set_crs(4326)

#then you can match the ID from gnv_latlon to 
gainsville_df <- fread("new_dataframe_data.csv", drop = c("Latitude","Longitude", "Census Code"))

gainsville_df <-  merge(gnv_latlon, gainsville_df, by = "ID")

#remove latitude and longitude points that fall outside of the polygon
dplyr::mutate(gainsville_df, check = as.vector(sf::st_intersects(x = gnv_latlon, y = gnv_poly, sparse = FALSE))) -> outliers_before
sf::st_filter(x= outliers_before, y= gnv_poly, predicate= st_intersects) -> gainsville_df

#Took out my census api key because of a feed back from a SO member. Please add a comment
#if you would like my census key.

#I use this function from tidycensus to retrieve the country shapfiles. 

alachua <- tidycensus::get_acs(state = "FL", county = "Alachua",  geography = "tract", geometry = T, variables = "B01003_001")
gainsville_df$Geomtry <- NULL
gainsville_df$Geomtry <- alachua$geometry[match(as.character(gainsville_df$`Geo ID`), alachua$GEOID)]

#gets us the first graph with bounry
ggplot() + 
  geom_sf(data = gainsville_df,aes(geometry= Geomtry, fill= Population), alpha= 0.2) +
  coord_sf(crs = "+init=epsg:4326")+ 
  geom_sf(data= gnv_poly) #with alpha added, we get the transparent boundary

现在我想在不做任何未来手动操作的情况下获得第二张图像。
由此.....

对此，可能吗？

发现这个比较空间多边形并保留或删除 R 中的公共边界， 但这里的人只想从一个 shapefile 中删除边界。我试图操纵它。

编辑这是我在 SymbolixAU 方向之后尝试过的，但我的idx变量是来自1:7

fl <- sf::st_read("PATH\\GIS_cgbound\\cgbound.shp") %>%  sf::st_transform(crs = 4326)
gainsville_df$Geomtry <-  sf::st_as_sf(gainsville_df$Geomtry) %>%  sf::st_transform(crs= 4326)

#normal boundry plot
plot( fl[, "geometry"] )

# And we can make a boundary by selecting some of the goemetries and union-ing them
boundary <- fl[ gnv_poly$geometry %in% gainsville_df$Geomtry, ]
boundary <- sf::st_union( fl ) %>% sf::st_as_sf()

## So now 'boundary' represents the area you want to cut out of your total shapes

## So you can find the intersection by an appropriate method
## st_contains will tell you all the shapes from 'fl' contained within the boundary
idx <- sf::st_contains(x = boundary, y = fl)

#doesn't work, thus no way of knowing the overlaps
#plot( fl[ idx[[1]], "geometry" ] ) 

#several more plots which i can't make sense of
plot( fl[ st_intersection(gainsville_df$Geomtry, gnv_poly$geometry), ])
plot(gainsville_df$Geomtry) #this just plots tracts

Answer 1

我要使用 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"
  )