[R-sig-Geo] writeRaster to TIF changes raster maximum value
Guillaume Drolet
droletguillaume at gmail.com
Tue Sep 17 20:13:01 CEST 2013
Hi Robert,
Following my previous post, I tried with a smaller extent than the one
that was used previously:
e <- extent(457785, 516302, 5147085, 5217408)
For the rest, I ran the same code as before:
> rc <- crop(r, e)
> NAvalue(rc) <- 0
> t <- writeRaster(rc, y, overwrite = TRUE, datatype = "INT1U", NAflag = 0)
> t
class : RasterLayer
dimensions : 2344, 1951, 4573144 (nrow, ncol, ncell)
resolution : 30, 30 (x, y)
extent : 457785, 516315, 5147085, 5217405 (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_sub\L71015027_02720000730_B10.tif
names : L71015027_02720000730_B10
values : 0, 255 (min, max)
Don't know why but the minimum is now 0, while for the extent we used
before (which was larger) the minimum value was 1...
There are definitely some behaviour I don't quite get in all this...
Thanks,
Guillaume
2013/9/17 <Guillaume.Drolet at mrn.gouv.qc.ca>:
>
>
>
> -------- Message original --------
> Sujet: Re: [R-sig-Geo] writeRaster to TIF changes raster maximum value
> Date : Tue, 17 Sep 2013 01:25:53 -0400
> De : Robert J. Hijmans <r.hijmans at gmail.com>
> Pour : Drolet, Guillaume (DRF) <Guillaume.Drolet at mrn.gouv.qc.ca>
> Copie à : r-sig-geo <r-sig-geo at r-project.org>
>
>
>
> Re: [R-sig-Geo] writeRaster to TIF changes raster maximum value
>
> 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:
> > 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 :
> >> 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:
> >> > 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
> >> > --
> >> > Guillaume Drolet ing.f M.Sc.
> >> > Chercheur en télédétection et dynamique forestière
> >> > Direction de la recherche forestière
> >> > Ministère des ressources naturelles du Québec
> >> > 2700, rue Einstein, C.RC.335
> >> > Québec (Québec), Canada G1P 3W8
> >> > T : +1 418 643 7994 # 6727
> >> > F : +1 418 643 2165
> >> > guillaume.drolet at mrn.gouv.qc.ca
> >> > http://www.mrn.gouv.qc.ca/forets
> >> >
> >> > _______________________________________________
> >> > R-sig-Geo mailing list
> >> > R-sig-Geo at r-project.org
> >> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >>
> >
> > --
> > Guillaume Drolet ing.f M.Sc.
> > Chercheur en télédétection et dynamique forestière
> > Direction de la recherche forestière
> > Ministère des ressources naturelles du Québec
> > 2700, rue Einstein, C.RC.335
> > Québec (Québec), Canada G1P 3W8
> > T : +1 418 643 7994 # 6727
> > F : +1 418 643 2165
> > guillaume.drolet at mrn.gouv.qc.ca
> > http://www.mrn.gouv.qc.ca/forets
>
>
More information about the R-sig-Geo
mailing list