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

[Agustin Lobo]

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).

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