Skip to content

writeRaster to TIF changes raster maximum value

4 messages · Guillaume.Drolet at mrn.gouv.qc.ca, Robert J. Hijmans

Original
Guillaume.Drolet at mrn.gouv.qc.ca
# 
Dear list members,

I need your help trying to figure out why my raster's maximum value gets
changed in the resulting TIF file from a call to the writeRaster function.

I read here
(http://r-sig-geo.2731867.n2.nabble.com/writeRaster-to-ascii-file-asc-td7584093.html#none) 

that this problem could be due to my version of GDAL (1.11dev released
2013/04/13). This is the GDAL version that is in my system's path. Don't
know if this can help but when I load rgdal, I can read this:

Loaded GDAL runtime: GDAL 1.9.2, released 2012/10/08

Here are the details of my raster (myRaster) before calling writeRaster:

  > raster(myRaster)

class       : RasterLayer
dimensions  : 7261, 7961, 57804821  (nrow, ncol, ncell)
resolution  : 30, 30  (x, y)
extent      : 457785, 696615, 5147085, 5364915  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
+ellps=WGS84 +towgs84=0,0,0
data source : E:\temp\Fmask\Image\L71015027_02720000730_B10.tif
names       : L71015027_02720000730_B10
values      : 0, 255  (min, max)


I crop "myRaster" using the crop function with a bounding box and then I 
save the cropped raster as a TIF with:

  > writeRaster(rc, "myCroppedRaster.tif", datatype = "INT2S")
rgdal: version: 0.8-11, (SVN revision 479M)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 1.9.2, released 2012/10/08
Path to GDAL shared files:
C:/Users/Utilisateur/Documents/R/win-library/3.0/rgdal/gdal
GDAL does not use iconv for recoding strings.
Loaded PROJ.4 runtime: Rel. 4.7.1, 23 September 2009, [PJ_VERSION: 470]
Path to PROJ.4 shared files:
C:/Users/Utilisateur/Documents/R/win-library/3.0/rgdal/proj

class       : RasterLayer
dimensions  : 3628, 3994, 14490232  (nrow, ncol, ncell)
resolution  : 30, 30  (x, y)
extent      : 457785, 577605, 5147085, 5255925  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
+ellps=WGS84 +towgs84=0,0,0
data source : E:temp\myRaster.tif
names       : myRaster
values      : 0, 254  (min, max)


As shown in the output from writeRaster, the maximum value in the output 
TIF file is now 254. I double-checked and
there are pixels with value 255 in the cropped area of the original
raster so the maximum value shouldn't be 254 but 255.

Here's the output of my sessionInfo():

R version 3.0.1 (2013-05-16)
Platform: x86_64-w64-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=French_Canada.1252  LC_CTYPE=French_Canada.1252
LC_MONETARY=French_Canada.1252
[4] LC_NUMERIC=C                   LC_TIME=French_Canada.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] rgdal_0.8-11  raster_2.1-49 sp_1.0-13

loaded via a namespace (and not attached):
[1] grid_3.0.1      lattice_0.20-15 tools_3.0.1

I thank you in advance for any hint that could help me solve this
problem.

Guillaume
Robert J. Hijmans
# 
I suspect that the original file has Byte values and that 255 is used,
perhaps erroneously, as NA. Can you give me access to the file?
Thanks, Robert
On Mon, Sep 16, 2013 at 8:10 AM, <Guillaume.Drolet at mrn.gouv.qc.ca> wrote:
Guillaume.Drolet at mrn.gouv.qc.ca
# 
Dear Robert,

Here's a link to 2 images: the original raster (full) and the cropped
one (sub):

https://www.dropbox.com/sh/8zo8fo964kr6ttk/pH1JIpTZ1C

The original raster has Byte values, as you suspected. But I don't know
how to handle that so that it's not used as NA value. From what I read, 
NA values in the original raster should be 0, and minimum and maximum 
values should be 1 and 255, respectively.

Thanks for helping,

Guillaume

Le 2013-09-16 11:41, Robert J. Hijmans a ?crit :

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Robert J. Hijmans
    

    

      
      #
    

  


  

      
Guillaume,

This indeed happens because of NA value handling problems.  The native
raster format uses 255 as the NA value for Byte files. That probably
needs to change. You write such a file because you crop a large area
and a temporary file is written. There are work arounds, see below,
but i will ponder changes to the package that will makes these
unnecessary.

The most direct way to avoid this is to do
library(raster)
r = raster('L71015027_02720000730_B10_full.tif')
e =extent( 457785, 577605, 5147085, 5255925)
x = crop(r, e, 'abc.tif', datatype='INT2S', overwrite=TRUE)


Alternatively,

y = crop(r, e)
NAvalue(x) = 0
y <- writeRaster(y, 'abcd.tif', datatype='INT2S', overwrite=TRUE)

Robert

        
On Mon, Sep 16, 2013 at 1:58 PM,  <Guillaume.Drolet at mrn.gouv.qc.ca> wrote: