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

[Agustin Lobo]

Tested with the Nan, gdalinfo -stat does not recalculate as expected.
Also, there is a problem
with summary()

> adls48_31_21
class       : RasterStack
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
names       : antofagastasierrav2.1, antofagastasierrav2.2,
antofagastasierrav2.3
min values  : 0, 0, 0
max values  : 6396, 8295, 8216

> summary(adls48_31_21)
Error: object 'sumobj' not found

> writeRaster(adls48_31_21,filename="adls48_31_21",format="GTiff",dataType="INT2U",NAflag=0,overwrite=T)
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
data source : /media/Iomega_HDD/VOLCAN/HYPERION/RHyperion/adls48_31_21.tif
names       : adls48_31_21.1, adls48_31_21.2, adls48_31_21.3
min values  : 148, 646, 1110
max values  : 24392, 12864, 13619

> GDALinfo("adls48_31_21.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        adls48_31_21.tif
apparent band summary:
   GDType Bmin  Bmax Bmean Bsd hasNoDataValue NoDataValue
1 Float32  148 24392   NaN NaN           TRUE           0
2 Float32  646 12864   NaN NaN           TRUE           0
3 Float32 1110 13619   NaN NaN           TRUE           0
Metadata:
AREA_OR_POINT=Area

alobo at delia:/media/Iomega_HDD/VOLCAN/HYPERION/RHyperion$ gdalinfo
-stats adls48_31_21.tif
Driver: GTiff/GeoTIFF
Files: adls48_31_21.tif
Size is 1181, 3251
Coordinate System is:
PROJCS["WGS 84 / UTM zone 19S",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-69],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",10000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","32719"]]
Origin = (643485.012701101019047,7151114.995386035181582)
Pixel Size = (29.974597798475870,-29.990772070132266)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  643485.013, 7151114.995) ( 67d34' 9.69"W, 25d45' 3.52"S)
Lower Left  (  643485.013, 7053614.995) ( 67d33'30.84"W, 26d37'51.92"S)
Upper Right (  678885.013, 7151114.995) ( 67d12'59.47"W, 25d44'49.50"S)
Lower Right (  678885.013, 7053614.995) ( 67d12'11.04"W, 26d37'37.34"S)
Center      (  661185.013, 7102364.995) ( 67d23'12.98"W, 26d11'21.01"S)
Band 1 Block=1181x1 Type=Float32, ColorInterp=Gray
  Min=148.000 Max=24392.000
  Minimum=148.000, Maximum=24392.000, Mean=nan, StdDev=nan
  NoData Value=0
  Metadata:
    STATISTICS_MAXIMUM=24392
    STATISTICS_MEAN=nan
    STATISTICS_MINIMUM=148
    STATISTICS_STDDEV=nan
Band 2 Block=1181x1 Type=Float32, ColorInterp=Undefined
  Min=646.000 Max=12864.000
  Minimum=646.000, Maximum=12864.000, Mean=nan, StdDev=nan
  NoData Value=0
  Metadata:
    STATISTICS_MAXIMUM=12864
    STATISTICS_MEAN=nan
    STATISTICS_MINIMUM=646
    STATISTICS_STDDEV=nan
Band 3 Block=1181x1 Type=Float32, ColorInterp=Undefined
  Min=1110.000 Max=13619.000
  Minimum=1110.000, Maximum=13619.000, Mean=nan, StdDev=nan
  NoData Value=0
  Metadata:
    STATISTICS_MAXIMUM=13619
    STATISTICS_MEAN=nan
    STATISTICS_MINIMUM=1110
    STATISTICS_STDDEV=nan

Agus

On Mon, Aug 27, 2012 at 2:02 PM, Agustin Lobo <alobolistas at gmail.com> wrote:
> Robert,
>
> On Mon, Aug 20, 2012 at 7:21 AM, Robert J. Hijmans <r.hijmans at gmail.com> wrote:
>> 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 understand "this is a design choice", but I think it is a wrong one.
> Having mean and sd
> equal to 0 (or NaN) is not the same as having unknown mean and sd. The
> point is that
> gdal (not only QGIS) calculate the stats if this information is not
> provided. If this info is set to 0,
> gdal and QGIS cannot know that it has to be recalculated.
>
>> 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?
> The correct values do help a lot, actually solve my problem.
> Have not been able to test for the case with NaN, but I'm pessimistic.
> There is no way
> of just not writing the info? gdalinfo (and qgis) would automatically
> calculate it if missing:
> "-stats
> Read and display image statistics. Force computation if no statistics
> are stored in an image.
> "
> Thanks,
>
> Agus