(lat, lon) WKT 坐标不能用 st_transform 很好地重新投影
Posted
技术标签:
【中文标题】(lat, lon) WKT 坐标不能用 st_transform 很好地重新投影【英文标题】：(lat, lon) WKT coordinates do not reproject well with 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 是的，我后来意识到，这是因为我将坐标分成三行以提高可见性，如果你把它们都放在同一行，它就可以工作。
【参考方案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()
【讨论】：
以上是关于(lat, lon) WKT 坐标不能用 st_transform 很好地重新投影的主要内容，如果未能解决你的问题，请参考以下文章