How to crop, project and merge raster data from either side of the antimeridian in R terra?

I am cropping data from global latlon rasters (EPSG: 4326) and then want to project that data into a suitable local coordinate reference system. The area I am cropping spans the antimeridian, so I want to crop the global data either side of the antimeridian, project those two 'halves' and then merge them back together once I have projected them.

However, when I project the 2 halves, they end up with slightly different resolutions. I can resample to get the same resolutions, and then merge them, but the resulting merged data is missing some values across the antimeridian.

I'd like to know why the resolutions of the 2 halves are different, and how to merge them without a resulting data gap.

Reprex below illustrates the question

library(terra)
#> terra 1.7.55
library(geodata)

#equal area projection for Fiji from projectionwizard.org
fiji_crs <- "+proj=laea +lon_0=-181.8896484 +lat_0=-17.73775 +datum=WGS84 +units=m +no_defs"

#get temperature data for Fiji
fiji_temp <- bio_oracle(path=tempdir(), "Temperature", "Mean", benthic=FALSE, time="Present")

#split data into left and right hand side of the antimeridian (just taking the first raster from the stack of rasters for this example)
fiji_temp_lhs <- crop(fiji_temp[[1]], ext(176, 180, -21, -12))
fiji_temp_rhs <- crop(fiji_temp[[1]], ext(-180, -177, -21, -12))

#project the two halves - they now have slightly different resolution
fiji_temp_list <- lapply(list(fiji_temp_lhs, fiji_temp_rhs), function(x) project(x, fiji_crs))

fiji_temp_list
#> [[1]]
#> class       : SpatRaster 
#> dimensions  : 109, 48, 1  (nrow, ncol, nlyr)
#> resolution  : 9164.313, 9164.313  (x, y)
#> extent      : -230080.8, 209806.2, -364282.8, 634627.3  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=laea +lat_0=-17.73775 +lon_0=-181.8896484 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#> source(s)   : memory
#> name        : Present.Surface.Temperature.Mean 
#> min value   :                         25.97523 
#> max value   :                         29.11796 
#> 
#> [[2]]
#> class       : SpatRaster 
#> dimensions  : 109, 37, 1  (nrow, ncol, nlyr)
#> resolution  : 9190.435, 9190.435  (x, y)
#> extent      : 196542.6, 536588.7, -368078.4, 633679.1  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=laea +lat_0=-17.73775 +lon_0=-181.8896484 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
#> source(s)   : memory
#> name        : Present.Surface.Temperature.Mean 
#> min value   :                         26.02962 
#> max value   :                         29.06818

# resample the RHS raster to have the same resolution as the LHS
fiji_temp_list[[2]] <- resample(fiji_temp_list[[2]], rast(extent = ext(fiji_temp_list[[2]]), res = res(fiji_temp_list[[1]]), crs = crs(fiji_temp_list[[2]])))

#can now merge the resulting rasters but there are missing values where the two halves of the data have been merged
fiji_temp_list |>
  sprc() |>
  merge() |>
  plot()

^{Created on 2023-11-09 with reprex v2.0.2}

Notes: I am cropping in 2 halves rather than using a single polygon/ bounding box because there is less computer processing involved: a single polygon crop results in a raster that spans the entire equator. I also want to avoid projecting the uncropped global data as that would also require a lot of memory/ processing time.

Solution

Then you have to manage somehow the artifacts caused by projection and resampling, when NA values are assigned. My attempt is to vectorize the raster first (to get the land surface), then crop, project and merge (your part of code) and use focal() function to replace NA with mean values and finally reapply land surface back, like:

library(terra)

fiji_crs <- "+proj=laea +lon_0=-181.8896484 +lat_0=-17.73775 +datum=WGS84 +units=m +no_defs"

fiji_temp <- geodata::bio_oracle(path=tempdir(), "Temperature", "Mean", benthic=FALSE, time="Present")

fiji_temp |>
  terra::plot()

Let's substitute NA to -999 and other values to 1 and vectorize it:

g <- subst(fiji_temp, NA, -999)
g[g>-999] <- 1

p <- terra::as.polygons(g) |>
  project(fiji_crs)

Crop, project, resample and merge:

fiji_temp_lhs <- crop(fiji_temp[[1]], ext(176, 180, -21, -12))
fiji_temp_rhs <- crop(fiji_temp[[1]], ext(-180, -177, -21, -12))

fiji_temp_list <- lapply(list(fiji_temp_lhs, fiji_temp_rhs), function(x) project(x, fiji_crs))

fiji_temp_list[[2]] <- resample(fiji_temp_list[[2]], rast(extent = ext(fiji_temp_list[[2]]), res = res(fiji_temp_list[[1]]), crs = crs(fiji_temp_list[[2]])))

h <- fiji_temp_list |>
  sprc() |>
  merge()

Apply the focal() function only to NA cells:

h <- focal(h, w=9, fun=mean, na.policy="only", na.rm=T)

plot(h)

Finally rasterize the polygon(s) and apply it to raster

i <- rasterize(p[2], h)
plot(h*i)

^{Created on 2024-01-04 with reprex v2.0.2}