Project XYZ data in swiss coordinates to WGS84 and plot

Question

my "long-term objective" is to plot a free topographic dataset of Switzerland (http://www.toposhop.admin.ch/de/shop/products/height/dhm25200_1) and create an overlay with a shapefile containing swiss borders and data points representing meteorological stations in R.

Plotting the shapefile with borders and add stations as points works well, but what has failed in many tries now is projecting the topography dataset in swiss coordinates to WGS84, in order to be able to plot it together with the borders and stations in WGS84.

What seems the best solution of what I've tried over the last days:

# 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"

Despite the ordering: is this error a result of round-off errors in the coordinate transformation? Does R provide a better solution to project and plot the data in terrain.colors() and in a a way that shape and points can be added to the plot?

Why I ask that directly is: The dataset is also available as ESRI ASCII Grid, but as I've tried to plot it in swiss coordinates (for a frist overview) the default color was red to yellow. I've tried to plot it with the image() function for terrain.colors(), but then neither points nor shape could not be added.

Can anyone help?

Thanks in advance!!!

Answer 1

You could read your data with rasterFromXYZ and then project this rastrer with projectRaster() to the desired CRS.

Answer 2

You can also use the functions of the cwhmisc package.

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

Answer 3

Unfortunately, the package cwhmisc is not so well documented. That's why I came up with the following.

• If you only have a single value to convert, you could use the website https://www.swisstopo.admin.ch/de/karten-daten-online/calculation-services/navref.html
• If you have a few values, the Geodetic REST Web services (REFRAME Web API) might be you choice (see code below) and https://www.swisstopo.admin.ch/de/karten-daten-online/calculation-services/m2m.html
• If you have many values, you can convert the values with only a very small error with the code below, see PDF refsys_e.pdf; chapter 4.2 "Approximate formulas for the direct transformation of Swiss projection coordinates to ellipsoidal WGS84 coordinates, 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

Data

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

Code Snippet 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"))

Output

$easting
[1] "7.811298488115001"

$northing
[1] "46.70310278713399"

Code Snippet 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_)

Output

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