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

【问题标题】：(lat, lon) WKT coordinates do not reproject well with st_transform(lat, lon) WKT 坐标不能用 st_transform 很好地重新投影
【发布时间】：2021-12-14 16:16:03
【问题描述】：

我在导入具有 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) 如果不可能，我该如何更改坐标的顺序，考虑到会有数百行并且并非所有多点特征的长度都相同。

【问题讨论】：

将点加载为 sf，将坐标提取为常规数据框字段 st_drop_geometry，然后使用 st_as_sf() 再次加载为 sf 对象，但坐标顺序是否正确？例如，tst_sf %>% mutate(lat = unlist(map(tst_sf$geometry,1)), lon = unlist(map(tst_sf$geometry,2))) %>% st_drop_geometry() %>% st_as_sf(coords = c("lon", "lat"))
@Skaqqs 遗憾地不适用于多点！这只会把它变成一个单点（POINT特征），其中第一个坐标是多点中第二个点的纬度，第二个坐标是多点中第一个点的纬度，POINT (52.85631 52.86765)
对我来说，第一步，tst 到 tst_sf 不起作用：ogr: corrupt data
将多点转换为点，但保留分组变量。使用我上面的建议，然后将基于组的点转换为多点作为最后一步？
@D.J 是的，我后来意识到，这是因为我将坐标分成三行以提高可见性，如果你把它们都放在同一行，它就可以工作。

标签： r spatial sf proj wkt

【解决方案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()

结果

编辑

更新上述答案以管理具有多行的数据集（参见下面的 cmets）

请在下面找到以下表示：

您的数据

# 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 (v2.0.1) 于 2021-11-09 创建}

【讨论】：

这行得通，但我必须让它成为一个循环/函数才能按行运行，因为当数据集有几行时它不起作用，每一行都是多点！如果您想添加它，很高兴接受它作为答案，否则当我有时间整理代码时，我会自己使用循环添加答案！
好的。我看到您在示例中添加了一些额外的行。因此，我将尝试修改代码，使其适用于不同的行。干杯
@Virginia Morera Pujol，我以为您在问题中提供的数据中添加了 MULTIPOINT 行。我错了！您是否可以在问题结束时进行编辑，以引入带有一些额外行的数据集，以便我可以进行一些测试。理想情况下，我会尝试为您提供矢量化解决方案，而不是for 循环；这样会更有效率。
你好@Virginia Morera Pujol。非常感谢您更新数据集。因此，请在上述答案的末尾找到我的编辑。我最终选择了for 循环，因为project() 函数的矢量化非常复杂。也就是说，我能够通过删除包含删除包含geometry 的旧列并添加新的geometry 列的步骤来优化代码。为此，我直接在for 循环中更新tst_sf 对象的Trail 列。干杯。

【解决方案2】：

这是一种蛮力方法，虽然不是最有效的，但似乎很有效。

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") %>% st_set_crs(4326)

#get the latitude & longitude
coordinates <- unlist(tst_sf$Trail)
lat <- coordinates [1:(length(coordinates)/2)]
lon <- coordinates [(length(coordinates)/2+1):length(coordinates)]

#rearrange the columns
#convert back into MULTIPOINT and 
tst_sf$Trail <-sfheaders::sfc_multipoint( matrix(c(lon, lat), ncol=2, byrow=FALSE) )
#redefine the CRS
tst_sf <- tst_sf %>% 
   st_as_sf(wkt = "Trail") %>% 
   st_set_crs(4326)

library(rnaturalearth)
world <- ne_countries(country = 'ireland', scale = "medium", returnclass = "sf")
#st_crs(world)
tst_sf2 <- tst_sf %>% st_transform(st_crs(world))
ggplot(tst_sf2) + 
   geom_sf(data = world) +
   geom_sf(col = "red") + 
   theme_bw()

【讨论】：