通过R中的sf将经度和纬度序列转换为多边形
Posted
技术标签:
【中文标题】通过R中的sf将经度和纬度序列转换为多边形【英文标题】：Convert sequence of longitude and latitude to polygon via sf in R 【发布时间】：2018-07-01 05:25:57 【问题描述】：我有五个经纬度组成这样的形状。
df <- c(order=1:5,
lon=c(119.4,119.4,119.4,119.5,119.5),
lat=c(-5.192,-5.192,-5.187,-5.187,-5.191))
如何使用 sf 这样的包轻松将它们转换为 sf 多边形数据框？
## Simple feature collection with 1 feature and 0 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 119.4 ymin: -5.192 xmax: 119.5 ymax: -5.187
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## geometry
## 1 POLYGON ((119.4 ...
【问题讨论】：
【参考方案1】：我看到这个问题出现在搜索结果中，所以我想我会提供一种更灵活的方法，从一系列lat 和lon 坐标中创建sf 中的多边形。
st_as_sf 有一个参数coords，它将把给定的点作为数据框中的坐标列，并将这些列转换为sfPOINT 几何。然后，因为sf 与dplyr 配合得很好，我们可以将st_combine 点转换为MULTIPOINT 和st_cast 以转换为POLYGON。与st_polygon 的“手动”构造相比，这样做的好处是我们不必仔细考虑关闭环或传递给构造函数的嵌套列表的正确级别，如果我们有超过一组坐标中的一个多边形，我们可以使用group_by 一次创建所有多边形。
注意从技术上讲，您可以在summarise 中使用do_union=FALSE 来做到这一点，但我认为这种语法更清晰一些，更类似于普通的summarise。
df <- data.frame(
lon = c(119.4, 119.4, 119.4, 119.5, 119.5),
lat = c(-5.192, -5.192, -5.187, -5.187, -5.191)
)
library(tidyverse)
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, proj.4 4.9.3
polygon <- df %>%
st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
summarise(geometry = st_combine(geometry)) %>%
st_cast("POLYGON")
polygon
#> Simple feature collection with 1 feature and 0 fields
#> geometry type: POLYGON
#> dimension: XY
#> bbox: xmin: 119.4 ymin: -5.192 xmax: 119.5 ymax: -5.187
#> epsg (SRID): 4326
#> proj4string: +proj=longlat +datum=WGS84 +no_defs
#> geometry
#> 1 POLYGON ((119.4 -5.192, 119...
plot(polygon)
由reprex package (v0.2.0) 于 2018 年 10 月 5 日创建。
【讨论】：
当我提供经纬度列表时，我似乎无法group_by。例如。对于两个多边形，我提供如下 df：df <- tibble(AC_NO = c(100, 101), lon = list(c(76.4, 76.38, 76.4, 76.38), c(76.4, 76.38, 76.4, 76.38)), lat = list(c(9.58, 9.58, 9.56, 9.56), c(9.55, 9.55, 9.53, 9.53))。然后：df %>% group_by(AC_NO) %>% st_as_sf( 等。我得到一个错误，“list”不能被强制输入“double”。关于我做错了什么有什么想法吗？
Necro，但看起来你有列表列 - 取消嵌套可能会解决它。我意识到我在上面不准确地使用了“列表”这个词，现在改了【参考方案2】：
相当于@Yo B. answer，但使用sf
library(sf)
df <- data.frame(lon=c(119.4,119.4,119.4,119.5,119.5),
lat=c(-5.192,-5.192,-5.187,-5.187,-5.191))
# You need first to close your polygon
# (first and last points must be identical)
df <- rbind(df, df[1,])
poly <- st_sf(st_sfc(st_polygon(list(as.matrix(df)))), crs = 4326)
poly
## Simple feature collection with 1 feature and 0 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 119.4 ymin: -5.192 xmax: 119.5 ymax: -5.187
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## st_sfc.st_polygon.list.as.matrix.df....
## 1 POLYGON ((119.4 -5.192, 119...
编辑以回答 cmets 中的问题
请参阅main sf vignette 以获得对sf、sfc 和sfg 对象的清晰详细的解释，总结如下：
用于表示简单特征的三个类是：
sf，具有特征属性和特征几何的表（data.frame），其中包含 sfc，包含每个特征（记录）的几何形状的列表列，由组成 sfg，单个简单特征的特征几何。
st_sfc 函数仅构建几何列（这是一个多边形列表 - 这里只有一个多边形）。 sfc 中的“c”代表“列”。函数st_sf 构建一个完整的sf 对象（它也有一个data.frame 类），它是一个带有几何列的数据框。在给定的示例中，多边形没有附加数据（无属性）。您可以通过构建 data.frame 来附加数据：
poly <- st_sf(data.frame(landuse = "Forest",
size = 23 ,
st_sfc(st_polygon(list(as.matrix(df))))),
crs = 4326)
poly
## ## Simple feature collection with 1 feature and 2 fields
## geometry type: POLYGON
## dimension: XYZ
## bbox: xmin: 1 ymin: 119.4 xmax: 5 ymax: 119.5
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## landuse size geometry
## 1 Forest 23 POLYGON Z ((1 119.4 -5.192,...
然后您可以从空间对象中提取这些元素中的每一个并检查它们的类：
完整的 sf 对象：带有 sfc 几何列的 data.frame
class(poly)
## "sf" "data.frame"
提取为列表的第三列：sfc 对象
class(poly[[3]])
## "sfc_POLYGON" "sfc"
几何列的第一个元素：sfg 多边形对象
class(poly[[3]][[1]])
## "XY" "POLYGON" "sfg"
【讨论】：
没有st_sf()，我也能完成同样的事情。例如：poly <- st_sfc(st_polygon(list(as.matrix(df))), crs = 4326) 谁能解释一下有什么区别（我不能）？【参考方案3】：
library(sfheaders) 在 CRAN 上从 20191004 可以获取一个 data.frame 并将其转换为 sf 对象
library(sf)
library(sfheaders)
df <- data.frame(
lon = c(119.4, 119.4, 119.4, 119.5, 119.5),
lat = c(-5.192, -5.192, -5.187, -5.187, -5.191)
)
sfheaders::sf_polygon(
obj = df
)
## given only two columns of data are in df there's no need to specify lon & lat arguments
# Simple feature collection with 1 feature and 1 field
# geometry type: POLYGON
# dimension: XY
# bbox: xmin: 119.4 ymin: -5.192 xmax: 119.5 ymax: -5.187
# epsg (SRID): NA
# proj4string:
# id geometry
# 1 1 POLYGON ((119.4 -5.192, 119...
【讨论】：
如何将此列字符串转换为多边形？ “列表（C（-72.6689780388482，-72.6697951984133，-72.670538307013，-72.6696918495337，-72.6689780388482，-37.8155440677891，-37.8164165487739，-37.815944377963，-37.8151004670485，-37.8155440677891））”跨度> 【参考方案4】：我不知道“sf”，但以防万一您指的是“sp”，这里是 SPDF 的完整构造
df <- data.frame(lon=c(119.4,119.4,119.4,119.5,119.5),
lat=c(-5.192,-5.192,-5.187,-5.187,-5.191))
require(sp)
spdf <- SpatialPolygonsDataFrame(
SpatialPolygons(
Srl=list(
Polygons(srl=list(
Polygon(coords=df)
), ID=1)
)
),
data=data.frame(a=1)
)
plot(spdf)
【讨论】：
以上是关于通过R中的sf将经度和纬度序列转换为多边形的主要内容，如果未能解决你的问题，请参考以下文章