Projecting Rasters, shift to the north. Can i correct it?-CodePudding

I would like to project one raster onto another and while doing so I think the values are changing their position "to the north".

Is this an expected behavior?

I was hoping to create a longlat raster to use it for lookups and GeoJSON generation.

Strangely (or maybe expected, I don't know) resulting GeoJSON positions are shifted (what feels like 10km) to the north.

Do I have a logical mistake somewhere?

This is an example:

x <- raster(ncol=900, nrow=900)
x_proj <- " proj=stere  lat_0=90  lat_ts=90  lon_0=10  k=0.93301270189  x_0=0  y_0=0  a=6378137  b=6356752.3142451802  to_meter=1000  no_defs "
proj <- CRS(x_proj)
extent(x) <- extent(-523.4622, 376.5378, -4658.645, -3758.645)
projection(x) <- x_proj
x[seq(450,455),seq(1,900)]<-1

new_raster<-raster(ncols=900,nrows=900)
new_raster_crs<- " proj=longlat  datum=WGS84  zone=34  no_defs  ellps=WGS84"
new_raster_proj <- CRS(new_raster_crs)
extent(new_raster) <- extent(3.5889,14.6209, 47.0705, 54.7405)
projection(new_raster) <- new_raster_proj
new_raster<-projectRaster(x,new_raster,method = "bilinear")

Plot of raster x

Plot of Raster new_raster

Is there something I could to with source/dest raster to create a "true" longlat lookup / GeoJSON possibility?

Is there a mistake somewhere ?

Can i maybe change y_0=0 value to correct this?

If thats the case how can I get the exact value of shift?

Currently I only see the change visually.

CodePudding user response：

That is as expected. Map projections distort (at least one of) shape, size, distance, and direction. In this case, you observe a change in shape.

You do make a mistake here:

new_raster_crs <- " proj=longlat  datum=WGS84  zone=34  no_defs  ellps=WGS84"

"zone" is only relevant for the "UTM" coordinate reference system (and perhaps others), and if you define a datum, you should not also define an ellipsoid. So it should be

new_raster_crs <- " proj=longlat  datum=WGS84"

But it seems that the other parts you add are simply ignored.

Another mistake is that you still use raster, as it has been replaced by terra. With terra it goes like this:

library(terra)
x <- rast(ncol=900, nrow=900, ext=c(-523.4622, 376.5378, -4658.645, -3758.645),
    crs=" proj=stere  lat_0=90  lat_ts=90  lon_0=10  k=0.93301270189  x_0=0  y_0=0  a=6378137  b=6356752.3142451802  to_meter=1000  no_defs")
    
x[seq(450,455),seq(1,900)]<-1

y <- rast(ncols=900, nrows=900, ext= c(3.5889,14.6209, 47.0705, 54.7405), crs=" proj=longlat  datum=WGS84")

z <- project(x, y)
plot(z)