1

私の「長期的な目標」は、スイスの無料の地形データセット ( http://www.toposhop.admin.ch/de/shop/products/height/dhm25200_1 ) をプロットし、スイスの国境を含むシェープファイルでオーバーレイを作成し、 R の気象観測所を表すデータ ポイント。

シェープファイルを境界線でプロットし、ステーションをポイントとして追加することはうまくいきますが、多くの試行で失敗したのは、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 に答える