无法转换缺少crs的sfc对象,但st_crs显示对象确实有crs

anauzrmj  于 2023-04-18  发布在  其他
关注(0)|答案(2)|浏览(428)

我正在尝试使用坐标系绘制一些经度/纬度坐标。

# Reprex
# Data I received was prepared this way 
library(sf)
library(ggplot2)
library(dplyr)
test_df <- 
  data.frame(id = LETTERS[1:6],
             lon = c(-71.2280234, -112.0863874, -80.2432839, -83.9638426, -76.9199742, -77.097455),
             lat = c(42.4286907, 33.5080943, 26.6374761, 43.4304646, 40.2398118, 39.0023604)) %>%
  sf::st_as_sf(coords = c('lon','lat')) %>% 
  rename(geo_col = geometry) %>% 
  as_tibble()

无CRS绘图

奇怪的是,如果我不设置CRS,我可以很好地绘制它。

new_df <- 
  test_df %>%
  mutate(geometry = geo_col) %>% 
  sf::st_as_sf() 

ggplot(new_df) + 
  geom_sf() + 
  theme_void()

使用CRS绘图

但是如果我尝试设置CRS,ggplot 2似乎表明CRS丢失,即使它明显存在。

# Plotting after setting CRS
new_df <- 
  test_df %>%
  mutate(geometry = geo_col) %>% 
  sf::st_as_sf() %>% 
  sf::st_set_crs(value = 4326) 

ggplot(new_df) + 
  geom_sf() + 
  theme_void()
Error in st_transform.sfc(X[[i]], ...) : 
  cannot transform sfc object with missing crs

但它确实有一个CRS...

> st_crs(new_df)
Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

以前的答案似乎建议您必须在完成任何转换之前设置crs。
Coordinate transformation of sf sfc object does not seem to be working
但是如果我在使用st_set_crs()之后再尝试使用st_transform(),这没有什么区别。
我想知道是否有关于我创建一个名为geometry的列,然后在其上运行st_as_sf()的内容?

mbyulnm0

mbyulnm01#

当你使用sf::st_as_sf()调用一个坐标向量时,你需要告诉你的R会话如何解释这些数字。它们是指球面上的度,还是平面上的米?或者你是一个美国人,感觉爱国并使用美国测量英尺?
我强烈怀疑lat & lon是球面上的度数,应该在WGS84的棱镜内解释,但你需要告诉你的R会话。
以下面这段代码为例,注意我对sf::st_as_sf()的调用做了一些调整:

library(sf)
library(dplyr)

test_df <- data.frame(id = LETTERS[1:6],
                      lon = c(-71.2280234, -112.0863874, -80.2432839, -83.9638426, -76.9199742, -77.097455),
                      lat = c(42.4286907, 33.5080943, 26.6374761, 43.4304646, 40.2398118, 39.0023604)) %>%
  sf::st_as_sf(coords = c('lon','lat'), 
               # this part is important!
               crs = 4326) %>% 
  rename(geo_col = geometry)
  
# to double check a visual overview; looks legit...
mapview::mapview(test_df)

unftdfkk

unftdfkk2#

如果,看起来是这样的,我误解了这个问题-问题不在于正确地将数据导入R,而是修复导入时由于任何原因而损坏的数据,但我们被它卡住了,因为这就是世界的运作方式:
还考虑下面的代码;它所做的是获取(坏,坏!)test_df tibble并将其转换回正确的{sf}对象。
它再次使用sf::st_as_sf,但这次指向geo_col作为包含几何图形的列,并在导入时将crs设置为4326。

library(sf)
library(ggplot2)
library(dplyr)
# misformed tibble to be repaired
test_df <- 
  data.frame(id = LETTERS[1:6],
             lon = c(-71.2280234, -112.0863874, -80.2432839, -83.9638426, -76.9199742, -77.097455),
             lat = c(42.4286907, 33.5080943, 26.6374761, 43.4304646, 40.2398118, 39.0023604)) %>%
  sf::st_as_sf(coords = c('lon','lat')) %>% 
  rename(geo_col = geometry) %>% 
  as_tibble() # this strips the magical juice! bad, very bad...

# new dataframe rises from ashes of the old!
new_df <- test_df %>% 
  st_as_sf(sf_column_name = "geo_col", crs = 4326)

# a visual check
mapview::mapview(new_df)

编辑:你的new_df的问题似乎是使用了一个mutate调用来构造新的几何列。因此,你的数据框有两个几何列:一个称为geometry和原始的geo_col。一次只有一个是活动的,并从“未知”解析为“已知”CRS(注意标题中的小度数符号,显示单位)。

new_df
# Simple feature collection with 6 features and 1 field
# Active geometry column: geo_col
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: -112.0864 ymin: 26.63748 xmax: -71.22802 ymax: 43.43046
# Geodetic CRS:  WGS 84
# A tibble: 6 × 3
#   id                 geo_col             geometry
# * <chr>          <POINT [°]>              <POINT>
# 1 A     (-71.22802 42.42869) (-71.22802 42.42869)
# 2 B     (-112.0864 33.50809) (-112.0864 33.50809)
# 3 C     (-80.24328 26.63748) (-80.24328 26.63748)
# 4 D     (-83.96384 43.43046) (-83.96384 43.43046)
# 5 E     (-76.91997 40.23981) (-76.91997 40.23981)
# 6 F     (-77.09745 39.00236) (-77.09745 39.00236)

在不同的CRS中有两个几何列在技术上是可行的,但通常不可取的,因为各种有趣的事情都可能发生...
恕我直言,这可能被认为是ggplot 2的一部分上的一个bug,它真的应该知道得更好--即使它被要求绘制一个sf对象,并且有一个名为geometry的列,它不是按照{sf}规范的活动几何列,它应该被忽略用于绘制;参见:

attr(new_df, "sf_column")
# [1] "geo_col"

但恐怕这是一个回避比解决更容易的错误。

相关问题