【问题标题】:R and GIS: Convert multiple lists of GPS coords within sf object into polygonsR和GIS:将sf对象内的多个GPS坐标列表转换为多边形
【发布时间】:2020-11-25 11:49:58
【问题描述】:

我有一个看起来像这样的 sf 对象:

structure(list(Value = c(0.0224448168716316, 0.0232573334948221, 
0.0227533509931156, 0.0214526198514539, 0.0202076447577734, 0.0198987318707359
), geometry = structure(list(structure(c(0.866271, 51.288072), class = c("XY", 
"POINT", "sfg")), structure(c(0.866264, 51.288062), class = c("XY", 
"POINT", "sfg")), structure(c(0.866258, 51.288052), class = c("XY", 
"POINT", "sfg")), structure(c(0.866252, 51.288042), class = c("XY", 
"POINT", "sfg")), structure(c(0.866245, 51.288032), class = c("XY", 
"POINT", "sfg")), structure(c(0.866239, 51.288022), class = c("XY", 
"POINT", "sfg"))), class = c("sfc_POINT", "sfc"), precision = 0, bbox = structure(c(xmin = 0.866239, 
ymin = 51.288022, xmax = 0.866271, ymax = 51.288072), class = "bbox"), crs = structure(list(
    epsg = 4326L, proj4string = "+proj=longlat +datum=WGS84 +no_defs"), class = "crs"), n_empty = 0L), 
    lon = c(`1` = 0.866271, `2` = 0.866264, `3` = 0.866258, `4` = 0.866252, 
    `5` = 0.866245, `6` = 0.866239), lat = c(`1` = 51.288072, 
    `2` = 51.288062, `3` = 51.288052, `4` = 51.288042, `5` = 51.288032, 
    `6` = 51.288022), bearing = c(203.698481252083, 200.617925511915, 
    200.617929623639, 203.698495141445, 200.617937896147, 200.617942031876
    ), p1 = structure(c(0.86626754319351, 0.866260971412482, 
    0.866254971412566, 0.866248543193856, 0.866241971412719, 
    0.866235971412793, 51.2880670617056, 51.2880569523547, 51.2880469523549, 
    51.2880370617061, 51.2880269523551, 51.2880169523552), .Dim = c(6L, 
    2L), .Dimnames = list(NULL, c("lon", "lat"))), p2 = structure(c(0.866291169366694, 
    0.866285120745778, 0.866279120739972, 0.866272169349123, 
    0.866266120728323, 0.866260120722494, 51.2880605588926, 51.2880512550893, 
    51.2880412550884, 51.2880305588895, 51.2880212550864, 51.2880112550854
    ), .Dim = c(6L, 2L), .Dimnames = list(NULL, c("lon", "lat"
    ))), p3 = structure(c(0.866298082980183, 0.866291177921382, 
    0.866285177915419, 0.86627908296192, 0.866272177903464, 0.866266177897481, 
    51.2880704354811, 51.2880613503796, 51.2880513503784, 51.288040435477, 
    51.2880313503759, 51.2880213503747), .Dim = c(6L, 2L), .Dimnames = list(
        NULL, c("lon", "lat"))), p4 = structure(c(0.86625083061034, 
    0.866242879232234, 0.866236879238075, 0.866231830627904, 
    0.866223879249719, 0.866217879255536, 51.2880834410928, 51.2880727448954, 
    51.2880627448964, 51.2880534410959, 51.2880427448984, 51.2880327448994
    ), .Dim = c(6L, 2L), .Dimnames = list(NULL, c("lon", "lat"
    ))), p5 = structure(c(0.866243916996385, 0.866236822056449, 
    0.866230822062455, 0.866224917014635, 0.866217822074403, 
    0.866211822080377, 51.2880735645039, 51.2880626496048, 51.2880526496061, 
    51.288043564508, 51.2880326496086, 51.2880226496098), .Dim = c(6L, 
    2L), .Dimnames = list(NULL, c("lon", "lat"))), p6 = structure(c(0.866291169342742, 
    0.866285120723053, 0.866279120717266, 0.866272169325172, 
    0.866266120705609, 0.866260120699783, 51.2880605588731, 51.288051255069, 
    51.2880412550681, 51.28803055887, 51.2880212550662, 51.2880112550652
    ), .Dim = c(6L, 2L), .Dimnames = list(NULL, c("lon", "lat"
    )))), sf_column = "geometry", agr = structure(c(Value = NA_integer_, 
lon = NA_integer_, lat = NA_integer_, bearing = NA_integer_, 
p1 = NA_integer_, p2 = NA_integer_, p3 = NA_integer_, p4 = NA_integer_, 
p5 = NA_integer_, p6 = NA_integer_), .Label = c("constant", "aggregate", 
"identity"), class = "factor"), row.names = c(NA, 6L), class = c("sf", 
"data.frame"))

我试图让所有点 lon 和 lats(例如 p1.lon、p1.lat、p2.lon、p2.lat、p3.lon、p3.lat 等)形成一个多边形df1 内的每一行

我尝试过使用 st_as_sf

df1 %>% select(c("p1","p2","p3","p4","p5","p6")) %>% 
  st_as_sf(coords = c('lon', 'lat'))

df1 %>% select(c("p1","p2","p3","p4","p5","p6")) %>% 
  mutate(st_as_sf(coords = c('lon', 'lat')))

df1  %>%  select(c("p1","p2","p3","p4","p5","p6")) %>%
  as.data.frame() %>% 
  st_as_sf( coords = c("lon", "lat"), crs = 4326 ) 

还有很多尝试将 p.lon 和 p.lat 值转换为 st_polygon 列表的变体

我应该如何引用 df1 中的 p1、p2、p3 等值,以及如何显示 st_polygon 这些是用来创建多边形的值? 我需要以某种方式更改 df1 的结构吗?

真的很想得到一些帮助。

【问题讨论】:

    标签: r list gis polygon sf


    【解决方案1】:

    df1 的结构似乎有点奇怪,所以我会质疑它最初是如何制作的,是否可以简化。


    但是,使用您当前的结构,稍加调整(使用data.table,因为它对我来说更容易使用),我们可以将数据放入“长”形式,然后使用sfheaders从 data.frame 制作多边形

    
    ## Add an 'id' so we know which points go into which polygons
    df1$id <- 1:nrow( df1 )
    
    dt <- as.data.table( df1 )
    
    dt <- melt(dt, id.vars = c("id"), measure.vars = names(dt)[grepl("\\.l", names( dt ) )])
    
    ## Because youre columns are named 'p1., p2., p3.' we can use the number to put the points
    ## in the correct order
    dt[, point_sequence := gsub("[^0-9]","",variable)]
    
    ## and remove the p1. bits
    dt[, lonlat := gsub("p[0-9]\\.", "", variable)]
    
    setorder( dt, id, point_sequence )
    
    ## reshape again to have separate lon & lat columns
    ## I'm doing this using a 'join' in data.table
    dt <- dt[ lonlat == "lon", .(id, point_sequence, value) ][
      dt[ lonlat == "lat", .(id, point_sequence, value) ]
      , on = .(id, point_sequence)
      , nomatch = 0
    ]
    
    setnames( dt, c("value", "i.value"), c("lon", "lat"))
    
    sf_poly <- sfheaders::sf_polygon(
      obj = dt
      , x = "lon"
      , y = "lat"
      , polygon_id = "id"
    )
    
    sf::st_crs( sf_poly ) <- 4326
    
    sf_poly
    # Simple feature collection with 6 features and 1 field
    # geometry type:  POLYGON
    # dimension:      XY
    # bbox:           xmin: 0.8662118 ymin: 51.28801 xmax: 0.8662981 ymax: 51.28808
    # geographic CRS: WGS 84
    # id                       geometry
    # 1  1 POLYGON ((0.8662675 51.2880...
    # 2  2 POLYGON ((0.866261 51.28806...
    # 3  3 POLYGON ((0.866255 51.28805...
    # 4  4 POLYGON ((0.8662485 51.2880...
    # 5  5 POLYGON ((0.866242 51.28803...
    # 6  6 POLYGON ((0.866236 51.28802...
    
    ## If you want the 'Value' column put back on
    sf::st_geometry( df1 ) <- NULL 
    
    sf_poly <- merge(
      x = sf_poly
      , y = as.data.frame( df1[, c("id", "Value")] )
      , by = "id"
    )
    

    【讨论】:

      猜你喜欢
      • 2019-03-11
      • 2020-11-14
      • 2019-10-09
      • 1970-01-01
      • 2023-02-02
      • 2021-08-17
      • 2021-09-02
      • 1970-01-01
      • 2020-12-29
      相关资源
      最近更新 更多