我想找到某個邊界框內的所有河流。我的最終目標是按名稱對它們進行子集化，以便我可以選擇要繪制的河流，但首先我必須知道我范圍內的要素名稱！我更喜歡使用 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"))
但首先我需要知道邊界框中特征的名稱。
uj5u.com熱心網友回復：
看起來你很接近。我添加了將 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')
Created on 2022-03-30 by the reprex package (v2.0.1)
轉載請註明出處，本文鏈接：https://www.uj5u.com/shujuku/453819.html