Re-projecting to Mollweide causes issues

231 Views Asked by At

I read in a raster, and take the bounding box of it.

library("terra")
library("sf")
library("maps")

my_box <- vect(st_as_sfc(st_bbox(template_rast), crs = crs(template_rast, proj=TRUE)))
maps::map('world')
plot(my_box, add=T)

enter image description here

I now want to crop another raster by this area. But the other raster has a different projection. Before doing the crop, I check that my box aligns with the same area. But it doesn't.

ghsl <- rast("GHS_POP_E2015_GLOBE_R2019A_54009_1K_V1_0.tif")

my_box_transformed <- project(my_box, crs(ghsl, proj=T))

plot(ghsl, col = "red")
plot(my_box_transformed, add=T)

I get it's something to do with projections, but I'm struggling with what exactly. Any ideas please?

enter image description here

Here are the crs info of the input layers

crs(template_rast, proj=TRUE)
"+proj=longlat +datum=WGS84 +no_defs"

crs(ghsl, proj=T)
"+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
1

There are 1 best solutions below

1
On

Your example is not fully reproducible so I cannot say why the area ends up where it does, but here is an example that works

Example data

library("terra")
library(geodata)
w <- world(path=".")
r <- rast(ext(13, 180, 27, 90), res=1)
v <- as.polygons(r, ext=T)

Project both

wm <- project(w, "+proj=moll")
vm <- project(v, "+proj=moll")

plot(wm)
lines(vm, col="red")

enter image description here

What might be happening is that your my_box has a longitude > 180. That is fine, in principle, but then you may need to first split the polygon before projecting.

box = rotate(my_box, 180, split=TRUE)