Undesired consequences of value replacement in raster

Question

I'm working with raster files based on remote sensing, see her for zipped example file: ftp://ftp.glcf.umd.edu/glcf/GLOPEM/10_Day/011_020/1992_011_020_npp_latlon/1992_011_020_npp_latlon.tif.gz.

For further analyses I'd like to change the cells with value -999 (i.e. water) to NA. I've tried:

for (i in 1:length(nlist)){
  r=raster(nlist[i])
  calc( r , function(x) { x[ x == -999 ] <- NA; return(x) } ,overwrite=T, filename=paste(names(r),"NA_edited"), format="GTiff")
}

f2 <- function(x) {
  v <- getValues(x)
  v[v<0]<-NA
  x <- setValues(x, v)
  return(x)
}

And

for (i in 1:length(nlist)){
  r=raster(nlist[i])
  r=f2(r)
  writeRaster(r,filename=paste(names(r),"NA_edited"), format="GTiff",overwrite=T)
}

where nlist is a vector of file locations (i.e. list.files()).

They do get rid of -999 and create NA values, but they also alters the other values for some reason. It seems like all values become slightly higher, e.g. max value of one layer goes from 160 to 177.

Is there some interpolation going on, or is there some mojo going on that I have missed out on? I'm not able to reproduce this on example rasters, such as raster(system.file("external/test.grd", package="raster")). Which makes me think it might have something to do with file formatting?

Cheers!

Answer 1

Since you're using the raster package, try using the function NAvalue to specify the NA values in your rasters instead of manually replacing the -999 values with NA's.

According to ?NAvalue:

NAvalue returns the value that is used to write NA values to disk (in ’raster’ type files). If you set the NA value of a Raster* object, this value will be interpreted as NA when reading the values from a file. Values already in memory will not be affected.

If the NA value is smaller than zero, all values smaller or equal to that number will be set to NA.

An example from the help file:

r1 <- raster(system.file("external/rlogo.grd", package="raster"))
r2 <- r1
NAvalue(r2)
NAvalue(r2) <- 255
#plot(r1)
#plot(r2)

For your problem, you can load the raster (say, r <- raster(file)), then use the function NAvalue(r) <- -999, and it should specify all values in the raster r with values <= -999 as NA.