将瑞士坐标中的 XYZ 数据投影到 WGS84 并绘图答案

【问题标题】：Project XYZ data in swiss coordinates to WGS84 and plot将瑞士坐标中的 XYZ 数据投影到 WGS84 并绘图
【发布时间】：2013-12-03 17:27:46
【问题描述】：

我的“长期目标”是绘制一个瑞士的免费地形数据集 (http://www.toposhop.admin.ch/de/shop/products/height/dhm25200_1)，并使用包含瑞士边界和代表 R 中气象站的数据点的 shapefile 创建一个叠加层。

绘制带有边界的 shapefile 并将站点添加为点效果很好，但现在在许多尝试中失败的是将瑞士坐标中的地形数据集投影到 WGS84，以便能够将其与边界和站点一起绘制WGS84。

过去几天我尝试过的最佳解决方案是什么：

# read xyz data:
topo=read.table("DHM200.xyz", sep=" ")
CH.topo.x= as.vector(topo$V1)
CH.topo.y= as.vector(topo$V2)
library(rgdal)
coord.topo <- data.frame(lon=CH.topo.x, lat=CH.topo.y)
coordinates(coord.topo) <- c("lon", "lat")
proj4string(coord.topo) <- CRS("+proj=somerc +lat_0=46.95240555555556+
lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel+
towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs") # CH1903 / LV03 (EPSG:21781)
CRS.new <- CRS("+init=epsg:4326")  # WGS84
coord.topo.wgs84 <- spTransform(coord.topo, CRS.new)
# this should have transformed the coordinates properly into a SpatialPoints object

# I now try to replace the old swiss coordinates by the degrees lat/lon:
topo[,1:2]=coordinates(coord.topo.wgs84)
topo=topo[order(topo$V1, topo$V2),]

# but creating a raster (that can be plotted!?)
topo.raster= rasterFromXYZ(topo, res=c(0.002516,0.00184), crs="+init=epsg:4326",
digits=5)
# returns the following error (either with or without "res" input):
# " x cell sizes are not regular"

尽管排序：这个错误是坐标变换中舍入错误的结果吗？ R 是否提供了更好的解决方案来投影和绘制terrain.colors() 中的数据，并且可以将形状和点添加到绘图中？

我直接问这个问题的原因是：数据集也可以作为 ESRI ASCII Grid 使用，但是当我尝试在瑞士坐标中绘制它时（对于第一个概述），默认颜色是红色到黄色。我尝试使用image() 函数为terrain.colors() 绘制它，但是点和形状都无法添加。

谁能帮忙？

提前致谢！！！

【问题讨论】：

您能否使用rasterFromXYZ 读取您的数据，并使用projectRaster() 将该栅格投影到新的CRS？
天哪！非常感谢，最好用！！！有时只是被直率蒙蔽了双眼；-) 当然，您可以将其发布为“真正的答案”！

标签： r plot projection coordinate-transformation wgs84

【解决方案1】：

不幸的是，cwhmisc 包没有很好的文档记录。这就是为什么我想出了以下内容。

如果您只有一个值要转换，您可以使用网站https://www.swisstopo.admin.ch/de/karten-daten-online/calculation-services/navref.html
如果您有几个值，您可能会选择 Geodetic REST Web 服务 (REFRAME Web API)（参见下面的代码）和https://www.swisstopo.admin.ch/de/karten-daten-online/calculation-services/m2m.html
如果你有很多值，你可以用下面的代码转换值，只有很小的错误，见PDF refsys_e.pdf；第 4.2 章“将瑞士投影坐标直接转换为椭球 WGS84 坐标的近似公式，URL：https://www.swisstopo.admin.ch/content/swisstopo-internet/en/online/calculation-services/_jcr_content/contentPar/tabs/items/documents_publicatio/tabPar/downloadlist/downloadItems/20_1467104436749.download/refsys_e.pdf

数据

df_chx_chy <- data.frame(chx=c(628500, 675500), chy=c(172500, 271500))

代码片段 1

library(httr)
library(jsonlite)
i <- 1
fromJSON(content(GET(paste0("http://geodesy.geo.admin.ch/reframe/lv03towgs84?easting=",df_chx_chy$chx[i], "&northing=",df_chx_chy$chy[i], "&format=json")), "text"))

输出

$easting
[1] "7.811298488115001"

$northing
[1] "46.70310278713399"

代码片段 2

library(dplyr)
df_chx_chy %>% mutate(
    y_ = (chx - 600000)/1000000,
    x_ = (chy - 200000)/1000000,
    lon = (2.6779094 +
               4.728982 * y_ +
               0.791484 * y_ * x_ +
               0.1306 * y_ * x_^2 -
               0.0436 * y_^3) * 100 / 36,
    lat = (16.9023892 +
               3.238272 * x_ -
               0.270978 * y_^2 -
               0.002528 * x_^2 -
               0.0447 * y_^2 * x_ -
               0.0140 *x_^3) * 100 / 36) %>% select(-y_, -x_)

输出

     chx    chy      lon      lat
1 628500 172500 7.811297 46.70310
2 675500 271500 8.442366 47.58985

【讨论】：

【解决方案2】：

您可以使用rasterFromXYZ 读取您的数据，然后使用projectRaster() 将此光栅投影到所需的CRS。

【讨论】：

【解决方案3】：

您还可以使用 cwhmisc 包的功能。

library(cwhmisc)
LB2YX(long, lat) # convert from long/lat to Swiss coordinates
YX2LB(yToEast, xToNorth) # convert from Swiss coordinates to long/lat

【讨论】：