I have ~100 different 4-band satellite imagery .tif rasters that have floating-point 32 bit depth. I need to convert these to 16-bit unsigned in R while scaling pixel values (not just trashing high values), but I have no idea where to begin for even a single raster, let alone the whole batch. Any help would be much appreciated!
Edit with more info: I searched the raster package documentation for bit, pixel, and depth keywords without much luck. Judging by min/max pixel values information found under dataType(), I want to go from FLT4S (32-bit floating point) to INT2U (16 bit). I tried setting datatype with writeRaster() as shown in the example there, but the output was just a black and white image rather than normal satellite imagery.
image
class : RasterStack
dimensions : 4300, 8909, 38308700, 4 (nrow, ncol, ncell, nlayers)
resolution : 3, 3 (x, y)
extent : 691032, 717759, 57492, 70392 (xmin, xmax, ymin, ymax)
crs : +proj=utm +zone=17 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
names : image.1, image.2, image.3, image.4
min values : 7.244503e-02, 8.278998e-02, 2.286164e-05, 8.571137e-02
max values : 0.5347134, 0.3522218, 0.4896736, 0.7308348
dataType(image) #[1] "FLT4S" "FLT4S" "FLT4S" "FLT4S"
image2 = writeRaster(image, 'new.tif', datatype='INT2U', overwrite=TRUE, format="GTiff")
dataType(image2) #[1] "INT2U"
image2
class : RasterBrick
dimensions : 4300, 8909, 38308700, 4 (nrow, ncol, ncell, nlayers)
resolution : 3, 3 (x, y)
extent : 691032, 717759, 57492, 70392 (xmin, xmax, ymin, ymax)
crs : +proj=utm +zone=17 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
source : new.tif
names : new.1, new.2, new.3, new.4
min values : 0, 0, 0, 0
max values : 1, 0, 0, 1
In your example you have real values between 0 and 1. By writing to an integer type, these are all truncated to 0. If you want to write to INT2U (unsigned byte) you could first scale the values to be between 0 and 255.
Example data
library(raster)
b <- brick(system.file("external/rlogo.grd", package="raster"))
image <- clamp(b/255, 0, 0.6)
#class : RasterBrick
#dimensions : 77, 101, 7777, 3 (nrow, ncol, ncell, nlayers)
#resolution : 1, 1 (x, y)
#extent : 0, 101, 0, 77 (xmin, xmax, ymin, ymax)
#crs : +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
#source : memory
#names : red, green, blue
#min values : 0, 0, 0
#max values : 0.6, 0.6, 0.6
What you do (truncation)
fname <- paste0(tempfile(), ".tif")
x <- writeRaster(image, fname, datatype='INT2U', overwrite=TRUE)
x
#min values : 0, 0, 0
#max values : 1, 1, 1
But note the difference when you use rounding
fname <- paste0(tempfile(), ".tif")
y <- round(image)
z <- writeRaster(y, fname, datatype='INT2U', overwrite=TRUE)
s <- stack(x[[1]], z[[1]])
plot(s)
Now with some scaling
maxv <- 65535
r <- round(image * maxv)
fname <- paste0(tempfile(), ".tif")
s <- writeRaster(r, fname, datatype='INT2U', overwrite=TRUE)
#s
#min values : 0, 0, 0
#max values : 39321, 39321, 39321
With your data you would get max values of
round(maxv * c(0.5347134, 0.3522218, 0.4896736, 0.7308348 ))
#[1] 35042 23083 32091 47895
You can also choose to set the max values of all layers to maxv, to keep more variation (but make the values no longer comparable with other data) --- more relevant if you were using a smaller range such as 0-255.
ss <- round(maxv * image / maxValue(image))
#names : red, green, blue
#min values : 0, 0, 0
#max values : 65535, 65535, 65535
the above works because the lowest values are zero; if you have negative values you would do
ss <- image - minValue(image)
ss <- round(maxv * ss / maxValue(ss))
In other cases you might want to use clamp. So you need to decide how to scale
What I showed is linear scaling. There are other ways. For example, there is the is also the scale method, to improve the statistical distribution of numbers. That may be relevant; but it depends on what your goals are.
Note. You do not say why you do this, and that is OK. But if it is to save hard disk space, you can use compression instead (see ?writeRaster)