通过边界框或多边形对空间特征进行子集化

Question

我想找到某个边界框内的所有河流。我的最终目标是按名称对它们进行子集化，这样我就可以选择绘制哪些河流，但首先我必须知道我范围内的要素名称！我更喜欢使用 ggplot/tidyverse 工具。

例如：

# Download river shapefile here: https://www.weather.gov/gis/Rivers

# Import river data as SF
st_read(dsn = 'rv16my07/', layer = 'rv16my07') %>%
  {. ->> my_rivers}

# Add a common CRS to the river dataset
st_crs(my_rivers) <-  CRS('+proj=longlat')

# Set x and y limits for the plot
ylims <- c(30.2, 31.4)
xlims <- c(-88.3, -87)

# Create sf for this bounding box
bounding.box <- st_as_sf(data.frame(long = xlims, lat = ylims), coords = c("lat", "long"), crs = CRS('+proj=longlat'))

# Try to find features that intersect
intersecting.rivers <- st_intersection(my_rivers, bounding.box)

然而，路口是空的。我在这里拧干什么？

在那个边界框内找到河流后，我想做这样的事情：

unqiue(intersecting.rivers$PNAME)

river.subsets <- subset(intersecting.rivers, PNAME %in% c("MOBILE R", "ALABAMA R"))

但首先我需要知道边界框内特征的名称。

Answer 1

看起来你很接近。我添加了将 xlim 和 ylim 中的两个对角点转换为 sf 对象并使用它们的边界框对河流进行子集化的缺失步骤。

library(sf)

# Download river shapefile here: https://www.weather.gov/gis/Rivers

# Import river data as SF
my_rivers <- st_read(dsn = '/home/x/Downloads/rivers/', layer = 'rv16my07') 
#> Reading layer `rv16my07' from data source `/home/x/Downloads/rivers' using driver `ESRI Shapefile'
#> Simple feature collection with 61122 features and 17 fields
#> Geometry type: MULTILINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -124.7068 ymin: 25.83636 xmax: -67.11324 ymax: 52.80121
#> Geodetic CRS:  NAD83

# The rivers data comes with a crs, so this step wasn't needed.
# Add a common CRS to the river dataset
#st_crs(my_rivers) <-  CRS('+proj=longlat')

# Set x and y limits for the plot, then make the points an sf object,
# set the crs as the same for my_rivers
ylims <- c(30.2, 31.4)
xlims <- c(-88.3, -87)
box_coords <- tibble(x = xlims, y = ylims) %>% 
  st_as_sf(coords = c("x", "y")) %>% 
  st_set_crs(st_crs(my_rivers))

#get the bounding box of the two x & y coordintates, make sfc
bounding_box <- st_bbox(box_coords) %>% st_as_sfc()


river_subset <- st_intersection(my_rivers, bounding_box)

head(river_subset)
#> Simple feature collection with 6 features and 17 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: -87.07318 ymin: 30.60249 xmax: -87 ymax: 30.83638
#> Geodetic CRS:  NAD83
#>      IHABBSRF_I          RR     HUC TYPE PMILE                  PNAME OWNAME
#> 8616       8616 03140104001 3140104    T   0.0           BLACKWATER R      0
#> 8617       8617 03140104002 3140104    R   0.5           BLACKWATER R      0
#> 8618       8618 03140104003 3140104    R   4.2           BLACKWATER R      0
#> 8630       8630 03140104015 3140104    R   7.8       BIG COLDWATER CR      0
#> 8631       8631 03140104016 3140104    R  17.3 BIG COLDWATER CR  E FK      0
#> 8634       8634 03140104019 3140104    R  17.3 BIG COLDWATER CR  W FK      0
#>           PNMCD OWNMCD       DSRR   DSHUC USDIR LEV J TERMID TRMBLV K
#> 8616 3140104001   <NA> 3140105007 3140105     R   1 0    205      1 0
#> 8617 3140104001   <NA> 3140104001 3140104     R   1 1    205      1 0
#> 8618 3140104001   <NA> 3140104002 3140104     R   1 1    205      1 0
#> 8630 3140104007   <NA> 3140104003 3140104     R   2 1    205      1 0
#> 8631 3140104008   <NA> 3140104015 3140104     R   2 2    205      1 0
#> 8634 3140104010   <NA> 3140104015 3140104     L   3 2    205      1 0
#>                            geometry
#> 8616 LINESTRING (-87.02298 30.60...
#> 8617 LINESTRING (-87.02928 30.60...
#> 8618 LINESTRING (-87.00626 30.64...
#> 8630 LINESTRING (-87 30.76254, -...
#> 8631 LINESTRING (-87.02458 30.78...
#> 8634 LINESTRING (-87.02458 30.78...

#Plot
ggplot() + 
  geom_sf(data = river_subset) + 
  geom_sf(data = bounding_box, fill = NA, color = 'red')

^{由 reprex package (v2.0.1)}

于 2022-03-30 创建

通过边界框或多边形对空间特征进行子集化

Subset Spatial Features by Bounding Box or Polygon

r

spatial

ggplot2

sf