Project XYZ data in swiss coordinates to WGS84 and plot

StackOverflow https://stackoverflow.com/questions/20357757

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!!!

Was it helpful?

Solution 2

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

OTHER TIPS

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

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

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
Licensed under: CC-BY-SA with attribution
Not affiliated with StackOverflow
scroll top