CODEWITHSUNDEEP
rplotprojectioncoordinate-transformationwgs84
user3017048
user3017048

Reputation: 2911

Project XYZ data in swiss coordinates to WGS84 and plot

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

Upvotes: 1

Views: 2602

Answers (3)

noslomo
noslomo

Reputation: 58

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

Upvotes: 1

takje
takje

Reputation: 2800

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

Upvotes: 0

EDi
EDi

Reputation: 13300

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

Upvotes: 0

Related Questions