6

我正在尝试将下面代码的输出（这是一个sfc 数据帧列表）转换为一个数据帧，其中列包含简单的特征集合- 每个观察一个多边形。

可重现示例的输出和以下预期输出：

# Reproducible example:  
library(tidyverse)  
library(sf)  
library(magrittr)   

# define radius for circle  
radius <- 40
r <- units::set_units(radius, units::as_units("nmile"), mode = "standard") 
%>% 
units::set_units(units::as_units("m"), mode = "standard")  

# Sample data: 
df <- data.frame(var = c("abc", "bcd", "cab", "dba"), 
lon = c(45,47,1, -109), 
lat = c(7, 10, 59, 30))

# Creating simple features with sf:
df <- df %>% st_as_sf(coords = c("lon", "lat"), dim = "XY")

# Applying Coordinate reference system WGS84:
df <- df %>% st_set_crs(4326)

# create function for finding UTM zones
utm_prj4 <- function(x) {
  coords <- cbind(x, st_coordinates(x))
  long <- coords$X
  lat <- coords$Y
  zone <- if(lat >= 56 && lat < 64 && long >= 3 && long < 12){x <- 32} else 
    if(
      lat >= 72 && lat < 84 && long >= 0 && long < 9) {x <- 31} else if(
        lat >= 72 && lat < 84 && long >= 9 && long < 21) {x <- 33} else if(
          lat >= 72 && lat < 84 && long >= 21 && long < 33) {x <- 35} else if(
            lat >= 72 && lat < 84 && long >= 33 && long < 42) {x <- 37} else{
              x <- (floor((long + 180)/6) %% 60) + 1
        }
  prj <- purrr::map2_chr(zone, lat, function(y, z){
    if (z >= 0){
      paste0("+proj=utm +zone=", y, " +datum=WGS84 +units=m +no_defs")
    } else{
      paste0("+proj=utm +zone=", y, " +south", " +datum=WGS84 +units=m 
         +no_defs")
    }})
  prj
}

# creates a list of data.frames, each with different crs
dfs <- map2(1:4, utm_prj4(df), function(x, y){
  st_transform(df[x,], y)
})
df <- map(dfs, ~ st_buffer(., r))

我试过 as.data.frame 但它不起作用。我试过 data.frame(t(as.data.frame(df)) 但它不起作用。

可重现示例的输出在左侧：

可重现示例的输出在左侧

4

1 回答 1

9

您需要确保您的sf对象具有相同的 crs，否则您无法将它们的几何图形组合成一个sfc（简单的特征列）。一旦我们将它们转换longlat为例如，我们就可以只使用rbind它们。

df_ll = map(df, ~ st_transform(., 4326))

df_sf = do.call(rbind, df_ll)
## or using Reduce
df_sf = Reduce(rbind, df_ll)
## or using purrr version of reduce
df_sf = reduce(df_ll, rbind)
于 2018-02-27T16:32:31.720 回答