[R-sig-Geo] wrong stats in GTiff files made with writeRaster() (through rgdal?)

Agustin Lobo alobolistas at gmail.com
Wed Aug 15 10:08:50 CEST 2012 

2012/8/14 Roger Bivand <Roger.Bivand at nhh.no>:
> On Tue, 14 Aug 2012, Agustin Lobo wrote:
>
>> I think this is a real need. The tif files currently generated by
>> writeRaster() cannot be easily modified
>> with gdal utilities, beacause gdalinfo can force the calculation of
>> the min/max values nut not of the
>> mean and sd unless they do not exist:
>>
>> from http://www.gdal.org/gdalinfo.html :
>>
.../...
>
>
> Is your problem that you want min, max, mean, sd in Qgis, but cannot
> generate them there?
Not exactly.  QGIS uses the stats for stretching if that info is
available in the file, and calculates it if not available.
But currently there is no way to force recalculation in case the stats
are wrong. Thus files created with writeRaster()
cannot be correctly displayed in QGIS.

> Have you tried using the setStatistics=TRUE argument in
> writeGDAL()?
No, that option is not documented in the help page.
I've tried now but the result is the same:
> writeRaster(adls48_31_21,filename="test",format="GTiff",dataType="INT2U",NAflag=0,overwrite=TRUE,setStatistics=TRUE )
class       : RasterBrick
dimensions  : 3251, 1181, 3839431, 3  (nrow, ncol, ncell, nlayers)
resolution  : 29.9746, 29.99077  (x, y)
extent      : 643485, 678885, 7053615, 7151115  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs
+ellps=WGS84 +towgs84=0,0,0
values      : /media/Iomega_HDD/VOLCAN/HYPERION/RHyperion/test.tif
min values  : 148, 646, 1110
max values  : 24392, 12864, 13619
layer names : test.1, test.2, test.3

> GDALinfo("test.tif")
rows        3251
columns     1181
bands       3
origin.x        643485
origin.y        7053615
res.x       29.9746
res.y       29.99077
ysign       -1
oblique.x   0
oblique.y   0
driver      GTiff
projection  +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs
file        test.tif
apparent band summary:
   GDType Bmin  Bmax Bmean Bsd hasNoDataValue NoDataValue
1 Float32  148 24392     0   0           TRUE           0
2 Float32  646 12864     0   0           TRUE           0
3 Float32 1110 13619     0   0           TRUE           0
Metadata:
AREA_OR_POINT=Area

> I don't think that writeRaster() passes it through, but does
> call RGDAL_SetStatistics internally in raster:::.stopGDALwriting(). If the
> statistics vector did not have 0, 0, for the mean and sd, but used any
> rgdal::writeGDAL flag, maybe the extra values could be output.
>
> If you read the code, you could have the pleasure of submitting patches!
which might be counteracted by your pleasure on patching my patches.
Let's wait until Robert  will be back.

Agus

> Roger
>
>
>
>>
>> By now I have to write the files again using gdal_translate -stats,
>> that forces stats calculation:
>>
>>
>> writeRaster(adls204_150_93,filename="adls204_150_93",format="GTiff",dataType="INT2U",NAflag=0,overwrite=T)
>> GDALinfo("adls204_150_93.tif") #wrong stats:
>> system("gdal_translate -stats adls204_150_93.tif test.tif")
>> system("gdal_translate -stats test.tif adls204_150_93.tif")
>> GDALinfo("adls204_150_93.tif")
>> system("rm test.tif")
>>
>> Agus
>>
>> 2012/8/14 Robert J. Hijmans <r.hijmans at gmail.com>:
>>>
>>> Agus,
>>> raster tracks the min and max values but does not track the mean and std;
>>> and most file formats cannot store them. The zeros are 'default' values;
>>> would it help to set them to another value that QGIS recognizes; perhaps
>>> to
>>> the NA value? I could add an option to compute and add the mean and std
>>> values.
>>> Robert
>>>
>>> On Fri, Aug 10, 2012 at 1:25 AM, Agustin Lobo <alobolistas at gmail.com>
>>> wrote:
>>>
>>>> I've just noted that files created with writeRaster() have wrong stats
>>>> (which is a problem later in Qgis, for example).
>>>> Not sure if writeRaster() uses rgdal, in that case the problem might
>>>> be more general:
>>>>
>>>>> delme <- subset(adlstest,c("B48","B31","B21"))
>>>>> summary(delme)
>>>>
>>>> Cells:  1089
>>>> NAs  :  0 0 0
>>>>
>>>>            1    2    3
>>>> Min.    2563 3491 3822
>>>> 1st Qu. 2700 3636 3971
>>>> Median  2764 3691 4009
>>>> Mean    2781 3699 4014
>>>> 3rd Qu. 2862 3765 4055
>>>> Max.    3629 3979 4258
>>>>
>>>>> GDALinfo("delme.tif")
>>>>
>>>> rows        33
>>>> columns     33
>>>> bands       3
>>>> origin.x        648490.8
>>>> origin.y        7058623
>>>> res.x       29.9746
>>>> res.y       29.99077
>>>> ysign       -1
>>>> oblique.x   0
>>>> oblique.y   0
>>>> driver      GTiff
>>>> projection  +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs
>>>> file        delme.tif
>>>> apparent band summary:
>>>>    GDType Bmin Bmax Bmean Bsd hasNoDataValue NoDataValue
>>>> 1 Float32 2563 3629     0   0           TRUE           0
>>>> 2 Float32 3491 3979     0   0           TRUE           0
>>>> 3 Float32 3822 4258     0   0           TRUE           0
>>>> Metadata:
>>>> AREA_OR_POINT=Area
>>>>
>>>> Note the 0 means and sd !
>>>>
>>>> adlstest here:
>>>> http://dl.dropbox.com/u/3180464/adlstest.rda
>>>>
>>>> Agus
>>>>
>>>>> sessionInfo()
>>>>
>>>> R version 2.15.1 (2012-06-22)
>>>> Platform: x86_64-pc-linux-gnu (64-bit)
>>>>
>>>> locale:
>>>>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
>>>> LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
>>>> LC_MONETARY=en_US.UTF-8
>>>>  [6] LC_MESSAGES=en_US.UTF-8    LC_PAPER=C                 LC_NAME=C
>>>>                LC_ADDRESS=C               LC_TELEPHONE=C
>>>> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
>>>>
>>>> attached base packages:
>>>> [1] stats     graphics  grDevices utils     datasets  methods   base
>>>>
>>>> other attached packages:
>>>> [1] rgdal_0.7-5   raster_2.0-08 sp_0.9-99
>>>>
>>>> loaded via a namespace (and not attached):
>>>> [1] grid_2.15.1     lattice_0.19-30 tools_2.15.1
>>>>
>>>> _______________________________________________
>>>> R-sig-Geo mailing list
>>>> R-sig-Geo at r-project.org
>>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>>>
>>>
>>>         [[alternative HTML version deleted]]
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo at r-project.org
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> --
> Roger Bivand
> Department of Economics, NHH Norwegian School of Economics,
> Helleveien 30, N-5045 Bergen, Norway.
> voice: +47 55 95 93 55; fax +47 55 95 95 43
> e-mail: Roger.Bivand at nhh.no
>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo

More information about the R-sig-Geo mailing list