调整R中sfc_POINT的值

ojsjcaue 于 2023-05-11 发布在其他

关注(0)|答案(2)|浏览(110)

对于美国，我可以得到一个州所有县的质心。然而，经过仔细检查，一些县的质心是不正确的。如何手动更正给定县的点几何图形（纬度和经度）？
弗吉尼亚州的例子

library(dplyr) # data wrangling
library(sf) # for getting centroids
library(tigris) # for getting county shapefile
library(cdlTools) # converting state fips to state name
library(ggplot2) # for mapping

### Extract centroids for each county in the USA

centroids <- counties(resolution = "20m") %>% 
  st_centroid()
centroids$STATE_NAME <- fips(centroids$STATEFP, to = 'Name')
va_cnty <- centroids[centroids$STATE_NAME %in% 'Virginia',]

### Make map of VA showing location of county centroids ----

state_shp <- states(cb = T)
state_shp$STATEFP <- as.numeric(state_shp$STATEFP)
state_shp <- state_shp[state_shp$STATEFP %in% c(1,3:14,16:56),]

cnty_shp <- counties(cb = TRUE)
cnty_shp$STATEFP <- as.numeric(cnty_shp$STATEFP)
lower48_shp <- cnty_shp[cnty_shp$STATEFP %in% c(1,3:14,16:56),]
va_shp <- lower48_shp[lower48_shp$STUSPS == 'VA',]

ggplot() +
  geom_sf(data = state_shp,
          mapping = aes(geometry = geometry),
          color = 'black',
          fill = 'gray70') +
  geom_sf(data = va_shp,
          mapping = aes(geometry = geometry),
          color = 'black',
          fill = 'lightyellow') +
  geom_sf(data = va_cnty,
          mapping = aes(geometry = geometry),
          color = "black",
          fill = 'red',
          shape = 21,
          size = 3) +
  coord_sf(xlim = c(-83.5, -75), ylim = c(36.4, 39.5), expand = TRUE) +
  theme_bw() +
  theme(text = element_text(size = 16),
        axis.text.x = element_text(size = 14, color = "black"),
        axis.text.y = element_text(size = 14, color = "black"),
        panel.grid.major = element_blank(),
        panel.background = element_rect(fill = "lightblue"))

弗吉尼亚州各县的中心经纬度Mapx1c 0d1x
如何更改va_cnty的几何列中的纬度和经度？正确的纬度和经度是：37.772236，-75.660827
两次不成功的尝试

va_cnty[va_cnty$NAME %in% 'Accomack',7] <- st_point(c(-75.660827,37.772236))

# Error in `[<-.data.frame`(`*tmp*`, va_cnty$NAME %in% "Accomack", 7, value = c(-75.660827,: replacement has 2 rows, data has 1

sfc = st_sfc(st_point(c(-75.660827,37.772236)), crs = 'NAD83')
va_cnty[va_cnty$NAME %in% 'Accomack',7] <- sfc

# Warning message:
# In `[<-.data.frame`(`*tmp*`, va_cnty$NAME %in% "Accomack", 7, value = # list( : replacement element 1 has 2 rows to replace 1 rows

来源：https://stackoverflow.com/questions/76204389/adjust-value-of-sfc-point-in-r

2条答案

按热度按时间

9lowa7mx1#

事实证明，您正在为弗吉尼亚州的县使用两个不同的shapefile。一个（counties(resolution = "20m")）返回一个覆盖整个县（包括海洋）的多边形。另一个（cnty_shp <- counties(cb = TRUE)）只返回县的陆地，至少在海中岛屿的例子中是这样。这导致非常不同的质心。
在前两个图中，请注意东部多边形的差异。第一个是陆地和海洋。在第二种方法中，只有陆地是由（多）多边形描述的
首先是counties(resolution = '20m')

virginia_counties <- counties(state = 'VA', resolution = '20m') 
virginia_centroids <- st_centroid(virginia_counties)

ggplot() + 
  geom_sf(data = virginia_counties) + 
  geom_sf(data = virginia_centroids, color = 'red')

所有的质心看起来都是正确的，因为它们在多边形中，但是其中一些多边形覆盖了非陆地区域。

第二个counties(cb = TRUE)

virginia_counties_2 <- counties(cb = TRUE, state = 'VA')
virginia_centroids_2 <- st_centroid(virginia_counties_2)

ggplot() + 
  geom_sf(data = virginia_counties_2) + 
  geom_sf(data = virginia_centroids_2, color = 'red')

两个质心在一起：

红色质心在具有等于（或几乎等于）围绕陆地的船体的陆地多边形的县中被重绘。

赞(0）回复(0）举报 2023-05-11

inn6fuwd2#

sf* 中的几何列是一个列表列，所以当处理单个元素时，你必须像处理列表一样处理它：

library(sf)
# sf example dataset
nc <- st_read(system.file("shape/nc.shp", package="sf")) |>
  st_centroid() 
nc[1:5,"NAME"]
#> Simple feature collection with 5 features and 1 field
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -81.49823 ymin: 36.40714 xmax: -76.02719 ymax: 36.49111
#> Geodetic CRS:  NAD27
#>          NAME                   geometry
#> 1        Ashe  POINT (-81.49823 36.4314)
#> 2   Alleghany POINT (-81.12513 36.49111)
#> 3       Surry POINT (-80.68573 36.41252)
#> 4   Currituck POINT (-76.02719 36.40714)
#> 5 Northampton POINT (-77.41046 36.42236)

# logical index 
nc$geometry[nc$NAME == "Ashe"] <- st_point(c(111,111))
# numeric index
nc$geometry[[2]] <- st_point(c(222,222))

# st_geometry method instead of $
st_geometry(nc)[[3]] <- st_point(c(333,333))

# replace just lon
nc$geometry[[4]][1] <- 123
# replace just lat
nc$geometry[[4]][2] <- 321

结果：

nc[1:5,"NAME"]
#> Simple feature collection with 5 features and 1 field
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -77.41046 ymin: 36.42236 xmax: 333 ymax: 333
#> Geodetic CRS:  NAD27
#>          NAME                   geometry
#> 1        Ashe            POINT (111 111)
#> 2   Alleghany            POINT (222 222)
#> 3       Surry            POINT (333 333)
#> 4   Currituck            POINT (123 321)
#> 5 Northampton POINT (-77.41046 36.42236)

创建于2023-05-09，使用reprex v2.0.2

赞(0）回复(0）举报 2023-05-11

我来回答

调整R中sfc_POINT的值

2条答案

相关问题

热门标签

最新问答