R projectRaster RGB Preservation – Preserving RGB Values When Using projectRaster in R

coordinate systemrraster

I'm using R and trying to change the CRS of a GeoTIFF (loaded as RasterStack) with RGB bands to a different CRS. To do this, I'm using projectRaster(). The issue I have is that it seems that the RGB values are being transformed too. Depending on the source and target CRS the RGB values can even become negative.

Here's an example:

library(raster)

# download sample GeoTIFF
download.file("https://github.com/EOxServer/autotest/raw/f8d9f4bde6686abbda09c711d4bf5239f5378aa9/autotest/data/meris/mosaic_MER_FRS_1P_RGB_reduced/mosaic_ENVISAT-MER_FRS_1PNPDE20060816_090929_000001972050_00222_23322_0058_RGB_reduced.tif", "sample.tif")

geotiff <- stack("sample.tif")

print(geotiff)

The output of that last command is the following (note the values of the three bands are in the range 0-255)

class      : RasterStack 
dimensions : 449, 541, 242909, 3  (nrow, ncol, ncell, nlayers)
resolution : 0.031355, 0.031355  (x, y)
extent     : 11.33175, 28.29481, 32.19025, 46.26864  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
names      : sample.1, sample.2, sample.3 
min values :        0,        0,        0 
max values :      255,      255,      255

Now, when I do the projection

geotiff_proj <- projectRaster(geotiff, crs = CRS("+proj=utm +zone=15 +datum=WGS84 +units=m +no_defs"))
print(geotiff_proj)

The output of the last command shows that the range of the RGB values has changed after the projection

class      : RasterBrick 
dimensions : 916, 744, 681504, 3  (nrow, ncol, ncell, nlayers)
resolution : 4140, 3220  (x, y)
extent     : 4813823, 7893983, 11461710, 14411230  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=15 +datum=WGS84 +units=m +no_defs 
source     : memory
names      : sample.1, sample.2, sample.3 
min values :        0,        0,        0 
max values : 252.7174, 253.5180, 253.4433

Is there a way to change the coordinate system of the raster without affecting the RGB values?

Best Answer

The transformation has to choose new RGB values for the cells in the transformed coordinate system grid. By default it uses a "bilinear" interpolation that smooths things a bit. I'm not sure how this can create values outside the range of the data (like your report of negative values, not reproduced here) unless possible some extrapolation might be doing it. Not sure.

The alternative method is "nearest neighbour". This is usually used for categorical data. Suppose 1=sea and 2=trees and 3=desert. A pixel that is half on trees and half in desert when recomputed via a transformation will get a value of 2.5 by interpolation, which is meaningless, so categorical data uses "nearest neighbour" which will always pick one value from the source grid cells for each destination grid cell.

If you do this for your RGB raster then you should get a raster which only has colours that are present in the source raster. This doesn't look much different to doing it by the default method. One of these is bilinear interpolation, the other is nearest neighbour. I'm not sure which is which:

enter image description here

geotiff_proj2 <- projectRaster(geotiff, 
    crs = CRS("+proj=utm +zone=15 +datum=WGS84 +units=m +no_defs"), 
     method="ngb")

You could also investigate the differences between these two rasters.

Related Question