Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

scale argument in gdal_translate() is wrongly built for multiband images #21

Open
aloboa opened this issue Nov 22, 2019 · 2 comments
Open

Comments

@aloboa
Copy link

aloboa commented Nov 22, 2019

#Goal: use gdal_transfom() to stretch floating valuesfrom 2% - 98%
quantiles to the 1- 255 range (to leave 0 for nodata) and write as
Byte (uint8)
#Test image:
m <- matrix(rnorm(n=10^4,m=0,sd=2),nrow=100)
m[1:10,] <- NA
imatest <- raster(m)
cellStats(imatest,range)
imatest <- brick(imatest,imatest*2,imatest+10)
extent(imatest) <- c(0,100,0,100)
projection(imatest) <- CRS("+init=epsg:25831")
imatest

class       : RasterBrick
dimensions  : 100, 100, 10000, 3  (nrow, ncol, ncell, nlayers)
resolution  : 1, 1  (x, y)
extent      : 0, 100, 0, 100  (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:25831 +proj=utm +zone=31 +ellps=GRS80
+towgs84=0,0,0,0,0,0,0 +units=m +no_defs
data source : in memory
names       :    layer.1,    layer.2,    layer.3
min values  :  -7.623476, -15.246953,   2.376524
max values  :   7.618672,  15.237344,  17.618672

writeRaster(imatest,filename="imatest",format="GTiff",NAflag=0,overwrite=TRUE)

#Calculate parameters for scale in the form:
#inmin1 inmax1 outmin1 outmax1 inmin2 inmax2 outmin2 outmax2 inmin3
inmax3 outmin3 outmax3

qt <- quantile(imatest,probs=c(0.02,0.98),na.rm=TRUE)
 qt

               2%       98%
layer.1 -4.060864  4.061431
layer.2 -8.121729  8.122863
layer.3  5.939136 14.061431

scal <- 1
for(i in 1:nlayers(imatest)){
    scal <- c(scal,c(qt[i,],c(1,255)))
}
scal <- scal[-1]
scal

        2%        98%                               2%        98%
                         2%
 -4.060864   4.061431   1.000000 255.000000  -8.121729   8.122863
1.000000 255.000000   5.939136
       98%
 14.061431   1.000000 255.000000

gdal_translate(src_dataset="imatest.tif",
               dst_dataset="imateststretchgdal.tif",
               ot="Byte",scale=scal)
# Warning message:
#   In system(cmd, intern = TRUE) :
#   running command '"C:\Program Files\QGIS
3.4\bin\gdal_translate.exe" -scale -4.03716844947439 4.10004776627608
1 255 -8.07433689894878 8.20009553255217 1 255 5.96283155052561
14.1000477662761 1 255 -ot "Byte" -of "GTiff" "imatest.tif"
"imateststretchgdal.tif"' had status 1
#
#Note scale is wrong, it should repeat "-scale" for each band
#-scale -4.03716844947439 4.10004776627608 1 255 -scale
-8.07433689894878 8.20009553255217 1 255 -scale 5.96283155052561
14.1000477662761 1 255
@jgrn307 jgrn307 changed the title scale argument in gdal_transform() is wrongly built for multiband images scale argument in gdal_translate() is wrongly built for multiband images Nov 22, 2019
@jgrn307
Copy link
Contributor

jgrn307 commented Nov 22, 2019

Can you see if the latest github version fixed the issue? And I assume this was for gdal_translate, not gdal_transform, right?

@aloboa
Copy link
Author

aloboa commented Dec 12, 2019

No, it's wrong. "-scale" has to be repeated once for each band, not for every parameter:


> gdal_translate(src_dataset="imatest.tif",
+                dst_dataset="imateststretchgdal.tif",
+                ot="Byte",scale=scal)
ERROR 6: Too many command options '5.82379595093665'
NULL
Warning message:
In system(cmd, intern = TRUE) :
  running command '"/usr/bin/gdal_translate" -scale -4.17620404906335 -scale 3.98042984386084 -scale 1 -scale 255 -scale -8.3524080981267 -scale 7.96085968772168 -scale 1 -scale 255 -scale 5.82379595093665 -scale 13.9804298438608 -scale 1 -scale 255 -ot "Byte" -of "GTiff" "imatest.tif" "imateststretchgdal.tif"' had status 1

it should be
gdal_translate imatest.tif imateststretchgdal2.tif -ot Byte -scale -4.17620404906335 3.98042984386084 1 255 -scale -8.3524080981267 7.96085968772168 1 255 -scale 5.82379595093665 13.9804298438608 1 255

Perhaps building the command is easier if scale is entered as matrix in R:

> mscal
               2%      98%      
layer.1 -4.176204  3.98043 1 255
layer.2 -8.352408  7.96086 1 255
layer.3  5.823796 13.98043 1 255

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants