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

[Roger Bivand]

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 :
>
> DESCRIPTION
>
> The gdalinfo program lists various information about a GDAL supported
> raster dataset.
>
> -mm
> Force computation of the actual min/max values for each band in the dataset.
> -stats
> Read and display image statistics. Force computation if no statistics
> are stored in an image.
>
> Perhaps a solution would be to let the files with "no statistics
> stored" instead of the current wrong
> default values. I think this would be ok for qgis also, as I think
> that qgis reads the mean and sd from
> the tif file if they are available (but with qgis you never know...).
> If possible, I would much prefer
> correct stat info than no info (but no info is better than wrong info).

Is your problem that you want min, max, mean, sd in Qgis, but cannot 
generate them there? Have you tried using the setStatistics=TRUE argument 
in writeGDAL()? 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!

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