将多边形形状文件导出为 XY 顶点对列表
Export a polygon shapefile as a list of XY vertex pairs
对于连续的多边形 shapefile,我需要一个格式如下的多边形顶点坐标文件作为遗留程序的输入:
polygon1name, AUTO
-120.408750 34.591250
-120.398313 34.591250
-120.396250 34.593313
-120.396250 34.593354
END
polygon2name, AUTO
-120.423354 34.641250
-120.423313 34.641250
-120.421250 35.643313
-120.421250 35.647521
END
从示例文件来看,遗留程序似乎需要逆时针绘图顺序中的对。下面以北卡罗来纳州各县为例。我希望在如何导出 XY 对和包含 , AUTO
和 END
部分方面得到帮助。
library(tidyverse) #for the %>% pipes and transmute()
library(sf) #for st_read()
library(rmapshaper) #for ms_simplify()
nc <- st_read(system.file("shape/nc.shp", package="sf")) %>%
transmute(NAME, geometry) %>% #keeps just the county column for simplicity
ms_simplify(keep = 0.01) #reduces the number of vertices for simplicity
plot(nc)
有什么想法吗?
谢谢。
不确定坐标的逆时针顺序,但这回答了您问题的输出格式部分。
#extract coordinates from sf
coord <- st_coordinates(nc) %>%
as.data.frame() %>%
group_by( L3 ) %>%
mutate(L4 = row_number() )
#extract data from sf
polygons <- st_drop_geometry(nc) %>%
mutate( NAME = as.character( NAME ) ) %>%
rownames_to_column( var = "id" ) %>%
mutate( id = as.numeric(id) ) %>%
#join coordinates
left_join( coord, by = c("id" = "L3") )
#split polygons-dataframe to list
l <- split( polygons, f = polygons$id )
#extract text needed from each polygon
result <- lapply( l, function(x) {
paste0 ( paste0( unique( x$NAME ), ", AUTO\n" ),
paste0( " ", x$X, " ", x$Y, collapse = "\n" ),
"\nEND" )
})
#unlist and write lines
writeLines( unlist(result) )
# Ashe, AUTO
# -81.4727554 36.2343559
# -81.7410736 36.3917847
# -81.6699982 36.5896492
# -81.3452988 36.5728645
# -81.2398911 36.3653641
# -81.4727554 36.2343559
# END
# Alleghany, AUTO
# -81.2398911 36.3653641
# -81.3452988 36.5728645
# -80.9034424 36.5652122
# -80.9563904 36.4037971
# -81.2398911 36.3653641
# END
# Surry, AUTO
# -80.4563446 36.2425575
# -80.874382 36.2338829
更新
对于逆时针部分:看一下sf::st_read()
的check_ring_dir
参数。
check_ring_dir
logical; if TRUE, polygon ring directions are checked
and if necessary corrected (when seen from above: exterior ring
counter clockwise, holes clockwise)
对于连续的多边形 shapefile,我需要一个格式如下的多边形顶点坐标文件作为遗留程序的输入:
polygon1name, AUTO
-120.408750 34.591250
-120.398313 34.591250
-120.396250 34.593313
-120.396250 34.593354
END
polygon2name, AUTO
-120.423354 34.641250
-120.423313 34.641250
-120.421250 35.643313
-120.421250 35.647521
END
从示例文件来看,遗留程序似乎需要逆时针绘图顺序中的对。下面以北卡罗来纳州各县为例。我希望在如何导出 XY 对和包含 , AUTO
和 END
部分方面得到帮助。
library(tidyverse) #for the %>% pipes and transmute()
library(sf) #for st_read()
library(rmapshaper) #for ms_simplify()
nc <- st_read(system.file("shape/nc.shp", package="sf")) %>%
transmute(NAME, geometry) %>% #keeps just the county column for simplicity
ms_simplify(keep = 0.01) #reduces the number of vertices for simplicity
plot(nc)
有什么想法吗?
谢谢。
不确定坐标的逆时针顺序,但这回答了您问题的输出格式部分。
#extract coordinates from sf
coord <- st_coordinates(nc) %>%
as.data.frame() %>%
group_by( L3 ) %>%
mutate(L4 = row_number() )
#extract data from sf
polygons <- st_drop_geometry(nc) %>%
mutate( NAME = as.character( NAME ) ) %>%
rownames_to_column( var = "id" ) %>%
mutate( id = as.numeric(id) ) %>%
#join coordinates
left_join( coord, by = c("id" = "L3") )
#split polygons-dataframe to list
l <- split( polygons, f = polygons$id )
#extract text needed from each polygon
result <- lapply( l, function(x) {
paste0 ( paste0( unique( x$NAME ), ", AUTO\n" ),
paste0( " ", x$X, " ", x$Y, collapse = "\n" ),
"\nEND" )
})
#unlist and write lines
writeLines( unlist(result) )
# Ashe, AUTO
# -81.4727554 36.2343559
# -81.7410736 36.3917847
# -81.6699982 36.5896492
# -81.3452988 36.5728645
# -81.2398911 36.3653641
# -81.4727554 36.2343559
# END
# Alleghany, AUTO
# -81.2398911 36.3653641
# -81.3452988 36.5728645
# -80.9034424 36.5652122
# -80.9563904 36.4037971
# -81.2398911 36.3653641
# END
# Surry, AUTO
# -80.4563446 36.2425575
# -80.874382 36.2338829
更新
对于逆时针部分:看一下sf::st_read()
的check_ring_dir
参数。
check_ring_dir
logical; if TRUE, polygon ring directions are checked and if necessary corrected (when seen from above: exterior ring counter clockwise, holes clockwise)