[R-sig-Geo] wrong stats in GTiff files made with writeRaster() (through rgdal?)
Roger Bivand
Roger.Bivand at nhh.no
Wed Aug 15 10:37:36 CEST 2012
On Wed, 15 Aug 2012, Agustin Lobo wrote:
> 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.
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.
> 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 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
>
--
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
More information about the R-sig-Geo
mailing list