Блог

У меня есть ~ 100 различных 4-полосных растров спутниковых изображений .tif, которые имеют 32-разрядную глубину с плавающей запятой. Мне нужно преобразовать их в 16-битный без знака в R при масштабировании значений пикселей (а не просто уничтожении высоких значений), но я понятия не имею, с чего начать даже для одного растра, не говоря уже о всей партии. Любая помощь будет высоко оценена!

Редактировать с дополнительной информацией: я искал в документации по пакету растров ключевые слова bit, pixel и depth без особой удачи. Судя по информации о минимальных / максимальных значениях пикселей, найденной в DataType(), я хочу перейти от FLT4S (32-разрядная плавающая точка) к INT2U (16 бит). Я попытался установить тип данных с помощью writeRaster (), как показано в примере, но на выходе было просто черно-белое изображение, а не обычные спутниковые снимки.

 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 
 

В вашем примере у вас есть реальные значения от 0 до 1. При записи в целочисленный тип все они усекаются до 0. Если вы хотите выполнить запись в INT2U (байт без знака), вы можете сначала масштабировать значения в диапазоне от 0 до 255.

Пример данных

 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 
 

Что вы делаете (усечение)

 fname <- paste0(tempfile(), ".tif")      
x <- writeRaster(image, fname, datatype='INT2U', overwrite=TRUE)
x
#min values : 0, 0, 0 
#max values : 1, 1, 1 
 

Но обратите внимание на разницу при использовании округления

 fname <- paste0(tempfile(), ".tif")      
y <- round(image)
z <- writeRaster(y, fname, datatype='INT2U', overwrite=TRUE)

s <- stack(x[[1]], z[[1]]) 
plot(s)
 

Теперь с некоторым масштабированием

 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 
 

С вашими данными вы получите максимальные значения

 round(maxv * c(0.5347134, 0.3522218, 0.4896736, 0.7308348 ))
#[1] 35042 23083 32091 47895
 

Вы также можете установить максимальные значения для всех слоев maxv , чтобы сохранить больше вариаций (но сделать значения несопоставимыми с другими данными) — более релевантными, если вы использовали меньший диапазон, например 0-255.

 ss <- round(maxv * image / maxValue(image))
#names      :   red, green,  blue 
#min values :     0,     0,     0 
#max values : 65535, 65535, 65535 
 

вышеуказанное работает, потому что наименьшие значения равны нулю; если у вас отрицательные значения, вы бы сделали

 ss <- image - minValue(image)
ss <- round(maxv * ss / maxValue(ss))
 

В других случаях вы можете использовать clamp . Итак, вам нужно решить, как масштабировать

То, что я показал, — это линейное масштабирование. Есть и другие способы. Например, существует также scale метод, позволяющий улучшить статистическое распределение чисел. Это может иметь значение; но это зависит от ваших целей.

Примечание. Вы не говорите, почему вы это делаете, и это нормально. Но если это необходимо для экономии места на жестком диске, вместо этого можно использовать сжатие (см. ?writeRaster )