您可以从包含顶点坐标的数据框中在 R 中创建多个多边形吗?
Can you create multiple polygons in R from a dataframe containing the vertices' coordinates?
我有一个包含 400 行的数据框,每行对应一个捕鱼区。每行包括所述捕鱼区域的顶点坐标,例如(NW = 西北角等):
> head(zones)
ID NW.X NW.Y NE.X NE.Y SE.X SE.Y SW.X SW.Y NW.X.2 Effort Season Method
1 1 -7,9961854 37.00222 -8,102379 37.05959 -8,1030245 37.01368 -7,9969111 36,9563306 -7,9961854 36 Winter Tremmel
2 2 -8,2268172 37.08076 -7,9773826 36.98974 -7,9778924 36.94389 -8,1248423 37,0183724 -8,2268172 12 Winter Gill
3 3 -8,102379 37.05959 -7,9961854 37.00222 -7,9990974 36.75608 -8,1036475 36,8128815 -8,102379 42 Winter Gill
4 4 -8,1024855 37.04711 -7,9963439 36.98973 -7,9994791 36.72262 -8,1036857 36,7794856 -8,1024855 42 Winter Tremmel
5 5 -8,3621785 37.02390 -7,9969701 36.93962 -7,9982745 36.83943 -8,3620979 36,9327666 -8,3621785 36 Winter Gill
6 6 -8,2580358 37.06539 -7,9963439 36.98973 -7,9974956 36.90627 -8,2580358 36,9822685 -8,2580358 36 Winter Gill
我想为每个单独的区域创建多边形,并将其提取为 shapefile 以在 ArcGIS 中使用。
我在这里看到了询问如何从坐标列表创建多边形的问题,但不是从以这种方式存储的坐标创建多个多边形的问题,我是 sf
包的初学者。
我知道通常情况下,您应该只有 1 列 X 和 1 列 Y,不像我这里的那样。我可以将我的数据框重新格式化为只有 2 列坐标,但是我该如何创建多个多边形?我想象有某种循环?
我认为这个问题的答案可能有帮助,但我不是很理解。 https://gis.stackexchange.com/questions/144728/add-polygons-to-spatialpolygons-via-loop-iteration-using-r
这是第一次尝试。不确定我是否不小心切换了 X 和 Y,所以请(目视)检查输出
library(data.table)
library(sf)
library(sfheaders)
library(tidyverse)
# Sample data -----
DT <- fread("ID NW.X NW.Y NE.X NE.Y SE.X SE.Y SW.X SW.Y NW.X.2 Effort Season Method
1 -7,9961854 37.00222 -8,102379 37.05959 -8,1030245 37.01368 -7,9969111 36,9563306 -7,9961854 36 Winter Tremmel
2 -8,2268172 37.08076 -7,9773826 36.98974 -7,9778924 36.94389 -8,1248423 37,0183724 -8,2268172 12 Winter Gill
3 -8,102379 37.05959 -7,9961854 37.00222 -7,9990974 36.75608 -8,1036475 36,8128815 -8,102379 42 Winter Gill
4 -8,1024855 37.04711 -7,9963439 36.98973 -7,9994791 36.72262 -8,1036857 36,7794856 -8,1024855 42 Winter Tremmel
5 -8,3621785 37.02390 -7,9969701 36.93962 -7,9982745 36.83943 -8,3620979 36,9327666 -8,3621785 36 Winter Gill
6 -8,2580358 37.06539 -7,9963439 36.98973 -7,9974956 36.90627 -8,2580358 36,9822685 -8,2580358 36 Winter Gill")
# Code -----
# Set decimal operator correct and convert all NW-like columns to numeric
cols <- grep("^(NW|NE|SE|SW)\.[XY]$", names(DT), value = TRUE)
DT[, (cols) := lapply(.SD, function(x) as.numeric(gsub(",", "\.", x))), .SDcols = cols]
#set to workable format df
my.poly <- setDF(DT) %>%
# Melt to long, beep XY paired
pivot_longer( cols = cols,
names_to = c("point", ".value"),
names_pattern = "(..)\.(.)" ) %>%
#!! check if X and Y are correct or should be switched!!
sfheaders::sf_polygon( x = "X", y = "Y", polygon_id = "ID" ) %>%
sf::st_set_crs(4326)
#visual incpection
mapview::mapview(my.poly)
我有一个包含 400 行的数据框,每行对应一个捕鱼区。每行包括所述捕鱼区域的顶点坐标,例如(NW = 西北角等):
> head(zones)
ID NW.X NW.Y NE.X NE.Y SE.X SE.Y SW.X SW.Y NW.X.2 Effort Season Method
1 1 -7,9961854 37.00222 -8,102379 37.05959 -8,1030245 37.01368 -7,9969111 36,9563306 -7,9961854 36 Winter Tremmel
2 2 -8,2268172 37.08076 -7,9773826 36.98974 -7,9778924 36.94389 -8,1248423 37,0183724 -8,2268172 12 Winter Gill
3 3 -8,102379 37.05959 -7,9961854 37.00222 -7,9990974 36.75608 -8,1036475 36,8128815 -8,102379 42 Winter Gill
4 4 -8,1024855 37.04711 -7,9963439 36.98973 -7,9994791 36.72262 -8,1036857 36,7794856 -8,1024855 42 Winter Tremmel
5 5 -8,3621785 37.02390 -7,9969701 36.93962 -7,9982745 36.83943 -8,3620979 36,9327666 -8,3621785 36 Winter Gill
6 6 -8,2580358 37.06539 -7,9963439 36.98973 -7,9974956 36.90627 -8,2580358 36,9822685 -8,2580358 36 Winter Gill
我想为每个单独的区域创建多边形,并将其提取为 shapefile 以在 ArcGIS 中使用。
我在这里看到了询问如何从坐标列表创建多边形的问题,但不是从以这种方式存储的坐标创建多个多边形的问题,我是 sf
包的初学者。
我知道通常情况下,您应该只有 1 列 X 和 1 列 Y,不像我这里的那样。我可以将我的数据框重新格式化为只有 2 列坐标,但是我该如何创建多个多边形?我想象有某种循环?
我认为这个问题的答案可能有帮助,但我不是很理解。 https://gis.stackexchange.com/questions/144728/add-polygons-to-spatialpolygons-via-loop-iteration-using-r
这是第一次尝试。不确定我是否不小心切换了 X 和 Y,所以请(目视)检查输出
library(data.table)
library(sf)
library(sfheaders)
library(tidyverse)
# Sample data -----
DT <- fread("ID NW.X NW.Y NE.X NE.Y SE.X SE.Y SW.X SW.Y NW.X.2 Effort Season Method
1 -7,9961854 37.00222 -8,102379 37.05959 -8,1030245 37.01368 -7,9969111 36,9563306 -7,9961854 36 Winter Tremmel
2 -8,2268172 37.08076 -7,9773826 36.98974 -7,9778924 36.94389 -8,1248423 37,0183724 -8,2268172 12 Winter Gill
3 -8,102379 37.05959 -7,9961854 37.00222 -7,9990974 36.75608 -8,1036475 36,8128815 -8,102379 42 Winter Gill
4 -8,1024855 37.04711 -7,9963439 36.98973 -7,9994791 36.72262 -8,1036857 36,7794856 -8,1024855 42 Winter Tremmel
5 -8,3621785 37.02390 -7,9969701 36.93962 -7,9982745 36.83943 -8,3620979 36,9327666 -8,3621785 36 Winter Gill
6 -8,2580358 37.06539 -7,9963439 36.98973 -7,9974956 36.90627 -8,2580358 36,9822685 -8,2580358 36 Winter Gill")
# Code -----
# Set decimal operator correct and convert all NW-like columns to numeric
cols <- grep("^(NW|NE|SE|SW)\.[XY]$", names(DT), value = TRUE)
DT[, (cols) := lapply(.SD, function(x) as.numeric(gsub(",", "\.", x))), .SDcols = cols]
#set to workable format df
my.poly <- setDF(DT) %>%
# Melt to long, beep XY paired
pivot_longer( cols = cols,
names_to = c("point", ".value"),
names_pattern = "(..)\.(.)" ) %>%
#!! check if X and Y are correct or should be switched!!
sfheaders::sf_polygon( x = "X", y = "Y", polygon_id = "ID" ) %>%
sf::st_set_crs(4326)
#visual incpection
mapview::mapview(my.poly)