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!!!
You could read your data with rasterFromXYZ and then project this rastrer with projectRaster() to the desired CRS.