3

I have a raster in UTM projection. The proj4 string is: +init=epsg:25833 +proj=utm +zone=33 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

I need to reproject the raster to DHDN / 3-degree Gauss-Kruger zone 4 (EPSG: 5678). So I did:

library('raster')

dtm <- read.table('333105624_dgm1.xyz', header = F, dec = '.')
dtm.raster <- rasterFromXYZ(dtm, res = c(1, 1), crs = crs('+init=epsg:25833'))
dtm.raster.rpj <- projectRaster(dtm.raster, crs = crs('+init=epsg:5678'), res = 1)

Unfortunatly, the reprojection doesn't align with the original:

The reprojected raster is shown above the original projected raster and as you can see it doesn't align with either the map or the original projected raster. What do I miss here?

2
  • 1
    Can you edit your quesiton and provide a reproducible example? That is create the input and expected output RasterLayers from scratch with raster(nrow=, ncol= ) etc. And perhaps include a minimal version of dtm, or at least show some of the values (you say projection is going bad, but perhaps dtm is not correct?). Can you also try this with terra rast(as.matrix(dtm), type="xyz") and project(dtm.raster) ? May 1 at 19:08
  • @RobertHijmans A link to the DTM in question (only about 100 mb in size) is included to the question (the word raster refers to the file). But anyway it works with the terra package. Thank you for that. It might be a workaround in my case but it doesn't deliver an answer why it doesn't work with the raster package.
    – Corbjn
    May 2 at 18:08
1

I only see a small difference between raster and terra. That is not unexpected as the extent is not exactly defined when projecting a raster based on a different CRS alone. It is better practice to provide a raster template to project to. Here is what I see:

options(rgdal_show_exportToProj4_warnings = "none")   
library(raster)
library(terra)

dtm <- read.table('333105624_dgm1.xyz', header = F, dec = '.')
dtm.raster <- rasterFromXYZ(dtm, crs = crs('+init=epsg:25833'))
dtm.raster.rpj <- projectRaster(dtm.raster, crs = crs('+init=epsg:5678'))


dtm.terra <- rast(as.matrix(dtm), type="xyz", crs ="epsg:25833")
dtm.terra.rpj <- project(dtm.terra, "epsg:5678")

plot(dtm.terra.rpj)
lines(extent(dtm.raster.rpj), col="red", lwd=4, xpd=T)
lines(ext(dtm.terra.rpj), col="blue", lwd=4, xpd=T)

enter image description here

e <- c(4522559, 4522565, 5624335, 5624341)
x <- crop(dtm.terra.rpj, ext(e))
y <- crop(dtm.raster.rpj, extent(e))

plot(y)
lines(as.polygons(x, dissolve=FALSE))

enter image description here

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.