r - (lat, lon) WKT 坐标不能用 st_transform 很好地重新投影

Question

我在导入具有 SRID 4326 的 wkt 多点特征的文件时遇到一些问题，其坐标是有序的(纬度、经度):

>st_crs(4326) 
Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
    ELLIPSOID["WGS 84",6378137,298.257223563,
        LENGTHUNIT["metre",1]]],
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(仅显示一行用于可重现的示例，但这将是数百行

tst <- data.frame(ID = rep("Test", 2), 
       SRIDTrail = rep(4326, 2), 
       Trail = c("MULTIPOINT (52.86 -8.00, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.89 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.97)", 
                 "MULTIPOINT (53.86 -7.00, 52.02 -6.98, 53.85 -7.80, 51.85 -8.98, 52.89 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.97)"))

tst_sf <- tst %>% 
  st_as_sf(wkt = "Trail") %>% 
  st_set_crs(4326)

现在，让我们从 naturalearth 包中下载世界地图，并检查它的 CRS:

library(rnaturalearth)
world <- ne_countries(scale = "medium", returnclass = "sf")
st_crs(world)

给出

Coordinate Reference System:
  User input: +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
  wkt:
BOUNDCRS[
    SOURCECRS[
        GEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]],
        CS[ellipsoidal,2],
            AXIS["longitude",east,
                ORDER[1],
                ANGLEUNIT["degree",0.0174532925199433,
                    ID["EPSG",9122]]],
            AXIS["latitude",north,
                ORDER[2],
                ANGLEUNIT["degree",0.0174532925199433,
                    ID["EPSG",9122]]]]],
TARGETCRS[
    GEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        CS[ellipsoidal,2],
            AXIS["latitude",north,
                ORDER[1],
                ANGLEUNIT["degree",0.0174532925199433]],
            AXIS["longitude",east,
                ORDER[2],
                ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]]],
ABRIDGEDTRANSFORMATION["Transformation from unknown to WGS84",
    METHOD["Geocentric translations (geog2D domain)",
        ID["EPSG",9603]],
    PARAMETER["X-axis translation",0,
        ID["EPSG",8605]],
    PARAMETER["Y-axis translation",0,
        ID["EPSG",8606]],
    PARAMETER["Z-axis translation",0,
        ID["EPSG",8607]]]]

表示第一个轴对应于经度，而不是纬度。所以，我尝试的第一件事(因为它会更快)是将我的数据转换为世界地图的相同投影，然后将它们都绘制出来:

tst_sf2 <- tst_sf %>% st_transform(st_crs(world))
ggplot(tst_sf2) + 
  geom_sf(data = world) +
  geom_sf(col = "red") + 
  theme_bw() 

这没有用，因为本应位于爱尔兰的点被绘制在印度洋中，该位置具有“交换”坐标，即纬度 -8，经度 53)。

让我们尝试反过来，变换世界地图，而不是 wtk。

world2 <- world %>% st_transform(st_crs(tst_sf))
ggplot(tst_sf) + 
  geom_sf(data = world2) +
  geom_sf(col = "red") + 
  theme_bw() 

这仍然不起作用:

所以，我的问题是:

(1) 是否有任何我可以使用的 EPSG 代码可以让 R 理解 WKT 文件中的坐标是根据预期交换的(我并不是要讨论应该是哪个顺序，只是为了修复它!)

(2) 如果这不可能，我该如何更改坐标的顺序，考虑到会有数百行并且并非所有多点要素的长度都相同。

Answer 1

请寻找另一个利用 sf_project() 函数的参数 authority_compliant = st_axis_order(FALSE/TRUE) 的解决方案。

• 您的数据

# The map
world <- ne_countries(scale = "medium", returnclass = "sf")

# Your point(s)
tst <- data.frame(ID = "Test",
                  SRIDTrail = 4326,
                  Trail = "MULTIPOINT (52.86 -8.00, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.89 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.97)")

tst_sf <- tst %>% 
  st_as_sf(wkt = "Trail", crs = 4326) # please, note that I lightened your code a little bit here.

• 代码

# set the new geometry (i.e. lon-lat instead of lat-long)
RightOrder <- sf_project(from = st_crs(tst_sf), 
                         to = st_crs(world), 
                         matrix(unlist(tst_sf$Trail), 
                                nrow = lapply(tst_sf$Trail, length)[[1]]/2, 
                                ncol = 2), 
                         authority_compliant = st_axis_order(TRUE)) %>% # the argument that allows to choose the order of the axes: lat-lon (FALSE) and lon-lat (TRUE) 
  as.data.frame() %>% 
  setNames(., c("lon", "lat")) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>% 
  st_cast("MULTIPOINT") %>% 
  st_union()

# drop the previous geometry and add the new one
tst_sf <- tst_sf %>%
  st_drop_geometry() %>%
  st_sf(.,RightOrder)

# visualize the result
ggplot(tst_sf) + 
  geom_sf(data = world) +
  geom_sf(col = "red") + 
  theme_bw() 

• 结果

---

编辑

更新上述答案以管理多行数据集(参见下面的评论)

请在下面找到以下代表:

• 您的数据

# The map
world <- ne_countries(scale = "medium", returnclass = "sf")

# Your point(s)
tst <- data.frame(ID = rep("Test", 2), 
                   SRIDTrail = rep(4326, 2), 
                   Trail = c("MULTIPOINT (52.86 -8.00, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.89 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.97)", 
                             "MULTIPOINT (53.86 -7.00, 52.02 -6.98, 53.85 -7.80, 51.85 -8.98, 52.89 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.98, 52.85 -7.97)"))
                  

tst_sf <- tst %>%  
  st_as_sf(wkt = "Trail", crs = 4326) # please, note that I lightened your code a little bit here.

• 代码

# set the new geometry (i.e. lon-lat instead of lat-long)
for (i in seq(tst_sf$Trail)){
  
  tst_sf$Trail[i] <- sf_project(from = st_crs(tst_sf), 
                                to = st_crs(world), 
                                matrix(unlist(tst_sf$Trail[i]), 
                                       nrow = lapply(tst_sf$Trail[i], length)[[1]]/2, 
                                       ncol = 2), 
                                authority_compliant = st_axis_order(TRUE)) %>% # the argument that allows to choose the order of the axes: lat-lon (FALSE) and lon-lat (TRUE) 
    as.data.frame() %>% 
    setNames(., c("lon", "lat")) %>% 
    st_as_sf(coords = c("lon", "lat"), crs = 4326) %>% 
    st_cast("MULTIPOINT") %>% 
    st_union()
}

# visualize the result
ggplot(tst_sf) + 
  geom_sf(data = world) +
  geom_sf(col = "red") + 
  theme_bw() 

• 结果

^{由 reprex package 创建于 2021-11-09 (v2.0.1)}