Search code examples
rgisr-rastermap-projections

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

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 x

Plot of Raster new_raster

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.

Solution

  • 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)