Best approach to minimize data loss when converting raster from projected to geographic coordinate system

0
votes

I have a worldwide raster file in PCS roughly at 1° resolution (extent is is 90 to -60 degrees latitude and 180 to -180 longitude) that I want to convert to GCS at 2° resolution but the current method that I use leads to data distortion. The origin raster contains categorical values but the same question would apply for continuous values. I wonder if there is a better way of doing this.

#create raster in PCS - CEA projection
x <- raster(ncol=360, nrow=142, xmn=-17367529, xmx=17367529, ymn=-6356742, ymx=7348382)
projection(x) <- "+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
values(x)=runif(51120, 1, 99)
# my target raster format
y <- raster(ncol=180, nrow=75, xmn=-180, xmx=180, ymn=-60, ymx=90)
projection(y) <- "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"
values(y)=runif(13500, 1, 99)
#downscale to match my target raster resolution
x2deg=raster::aggregate(x,fact=2)
#when I try to project I get the first error
x2deg.gcs = projectRaster(x2deg, y)

Error in if (maxy == miny) { : missing value where TRUE/FALSE needed In addition: Warning message: In rgdal::rawTransform(projfrom, projto, nrow(xy), xy[, 1], xy[, : 48 projected point(s) not finite

#So I slightly cut the extent and get the desired result
extent(x2deg) <- c(xmin= -17367529, xmax= 17367529, ymin= 0.999*(-6356742), ymax= .999*(7348382))
x2deg.gcs = projectRaster(x2deg, y)

Should I perhaps downscale ncol and nrows in raster x using a different factor? I know I will have to do some uneven interpolation to go from 142 rows to 75 but for example aggregate doesn't take fractions in the argument frac. Maybe there is another package that handles this better such as rgdal?

1
rgisrasterr-rasterrgdal

1 Answers

1
votes

This works better with terra

library(terra)
x <- rast(ncol=360, nrow=142, xmin=-17367529, xmax=17367529, ymin=-6356742, ymax=7348382)
crs(x) <- "+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs"
values(x) <- 1:ncell(x)

y <- rast(ncol=180, nrow=75, xmin=-180, xmax=180, ymin=-60, ymax=90)
crs(y) <- "+proj=longlat +datum=WGS84"

p <- project(x, y)