Agus,

As Roger suspected, this is a design choice; although possible, it is a
costly to compute the mean and sd when writing in chunks (the sd in
particular, as you need to have the mean to compute the sd) and most file
formats cannot store this anyway, and I was not aware of any application
using this.

I have changed raster's behavior (on R-Forge now) such that sd (denominator
is n-1) and mean are set to their correct values for Raster objects that
have the cell values in memory before writing the file, and to NaN (at
least that is what GDALinfo reports), instead of zero, for other cases
(when writing in chunks from another file). Does that help? That is, is the
NaN understood by QGIS?

Robert

On Wed, Aug 15, 2012 at 1:37 AM, Roger Bivand <Roger.Bivand@nhh.no> wrote:

> On Wed, 15 Aug 2012, Agustin Lobo wrote:
>
>  2012/8/14 Roger Bivand <Roger.Bivand@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<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.
>>
>
> Because - if you read the code - the mean and sd are set to 0.0 before
> writing. It isn't "wrong", it is a design choice not to incur the cost of
> calculating them by default, I think.
>
>
>
>>  Have you tried using the setStatistics=TRUE argument in
>>> writeGDAL()?
>>>
>> No, that option is not documented in the help page.
>>
>
> It has been documented since November 2010, for example bottom of p. 21 on
> http://cran.r-project.org/web/**packages/rgdal/rgdal.pdf<http://cran.r-project.org/web/packages/rgdal/rgdal.pdf>
> .
>
>
>
>  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 )
>>>
>>
> I did say:
>
>
>
>>  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!
>>>
>>
> However, if you convert to a SpatialGridDataFrame and use writeGDAL(),
> they are computed and displayed. I'm unsure from your case whether the
> missing values are correctly represented - I think all data were present -
> but the mean and sd are present.
>
>
>  which might be counteracted by your pleasure on patching my patches.
>> Let's wait until Robert  will be back.
>>
>
> Contributing code is a very important learning channel for the
> contributor, because the questioner/contributor is then in a much better
> position to understand their own questions, and possibly to contribute back
> something of value to others.
>
> Roger
>
>
>
>> 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@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@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<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@r-project.org
>>>>>> https://stat.ethz.ch/mailman/**listinfo/r-sig-geo<https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
>>>>>>
>>>>>>
>>>>>         [[alternative HTML version deleted]]
>>>>>
>>>>> ______________________________**_________________
>>>>> R-sig-Geo mailing list
>>>>> R-sig-Geo@r-project.org
>>>>> https://stat.ethz.ch/mailman/**listinfo/r-sig-geo<https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
>>>>>
>>>>
>>>>
>>>> ______________________________**_________________
>>>> R-sig-Geo mailing list
>>>> R-sig-Geo@r-project.org
>>>> https://stat.ethz.ch/mailman/**listinfo/r-sig-geo<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@nhh.no
>>>
>>>
>>> ______________________________**_________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo@r-project.org
>>> https://stat.ethz.ch/mailman/**listinfo/r-sig-geo<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@nhh.no
>
> ______________________________**_________________
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/**listinfo/r-sig-geo<https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
>

	[[alternative HTML version deleted]]

