找到点所在的状态,`sf` r
find the state where the points are located in, `sf` r
我有很多从卫星图像生成的点,我想找到这些点在哪里(状态)。我搜索了一些链接并得到了我应该使用 sf::st_intersects
的线索,但事实证明它不起作用。这是一个最小的例子:
library(sf)
library(ggplot2)
library(dplyr)
# Contiguous US state boundaries
usa = st_as_sf(maps::map("state", fill = TRUE, plot = FALSE))
# Simulate some random points
pts = data.frame(
x = c(-91.6, -74.3, -101.5),
y = c(36.1, 42.1, 25.3)
) %>%
st_as_sf(coords=c("x", "y"), crs = 4326)
ggplot() +
geom_sf(data = usa, fill = NA) +
geom_sf(data = pts,
shape = 21, size = 4, fill = "red") +
coord_sf(crs = st_crs(102003)) +
theme_minimal()
这是结果图:
我想要的是 pts data.frame 还有一行表示点的状态:
geometry state
1 POINT (-91.6 36.1) arkansas
2 POINT (-74.3 42.1) new york
3 POINT (-101.5 25.3) NA
我知道我应该显示 sf::st_transform
,但没有成功。理想情况下,我希望相交是可扩展的,因为我有超过 1,000,000,000 个点。
谢谢!
int = st_intersects(pts, usa)
sapply(int, function(x) if(length(x) == 0) NA else as.character(usa$ID[x]))
#[1] "arkansas" "new york" NA
您可以使用 st_join
.
a <- pts %>%
st_join(usa)
> a
Simple feature collection with 3 features and 1 field
geometry type: POINT
dimension: XY
bbox: xmin: -101.5 ymin: 25.3 xmax: -74.3 ymax: 42.1
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
ID geometry
1 arkansas POINT (-91.6 36.1)
2 new york POINT (-74.3 42.1)
3 <NA> POINT (-101.5 25.3)
我有很多从卫星图像生成的点,我想找到这些点在哪里(状态)。我搜索了一些链接并得到了我应该使用 sf::st_intersects
的线索,但事实证明它不起作用。这是一个最小的例子:
library(sf)
library(ggplot2)
library(dplyr)
# Contiguous US state boundaries
usa = st_as_sf(maps::map("state", fill = TRUE, plot = FALSE))
# Simulate some random points
pts = data.frame(
x = c(-91.6, -74.3, -101.5),
y = c(36.1, 42.1, 25.3)
) %>%
st_as_sf(coords=c("x", "y"), crs = 4326)
ggplot() +
geom_sf(data = usa, fill = NA) +
geom_sf(data = pts,
shape = 21, size = 4, fill = "red") +
coord_sf(crs = st_crs(102003)) +
theme_minimal()
这是结果图:
geometry state
1 POINT (-91.6 36.1) arkansas
2 POINT (-74.3 42.1) new york
3 POINT (-101.5 25.3) NA
我知道我应该显示 sf::st_transform
,但没有成功。理想情况下,我希望相交是可扩展的,因为我有超过 1,000,000,000 个点。
谢谢!
int = st_intersects(pts, usa)
sapply(int, function(x) if(length(x) == 0) NA else as.character(usa$ID[x]))
#[1] "arkansas" "new york" NA
您可以使用 st_join
.
a <- pts %>%
st_join(usa)
> a
Simple feature collection with 3 features and 1 field
geometry type: POINT
dimension: XY
bbox: xmin: -101.5 ymin: 25.3 xmax: -74.3 ymax: 42.1
epsg (SRID): 4326
proj4string: +proj=longlat +datum=WGS84 +no_defs
ID geometry
1 arkansas POINT (-91.6 36.1)
2 new york POINT (-74.3 42.1)
3 <NA> POINT (-101.5 25.3)