1

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

绘制带有边界的 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()for 函数绘制它terrain.colors()，但是点和形状都无法添加。

任何人都可以帮忙吗？

提前致谢！！！

4

3 回答 3

1

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

数据

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
于 2021-01-12T15:12:13.357 回答
0

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

于 2013-12-04T11:26:39.360 回答
0

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

library(cwhmisc)
LB2YX(long, lat) # convert from long/lat to Swiss coordinates
YX2LB(yToEast, xToNorth) # convert from Swiss coordinates to long/lat
于 2015-08-31T10:14:25.767 回答