Raster file statistics in R (stars and raster package) don't match gdalinfo

Question

I'm trying to read a gridded binary file in R using read_stars. Comparing the stats from stars package to the stats from gdalinfo, I could see differences in min/max values. Below is the sample code.

library(stars)
library(sf)
library(tidyverse)
library(fs)
library(raster)

url <- 'ftp://ftp.wpc.ncep.noaa.gov/pwpf/conus_2.5km/2.5kmpwpf_72hr/2.5kmprcntil_72hsnow_90pt_latestf072.grb'
file <- 'snow.grb'
download.file(url, destfile='snow.grb')

snow.stars <- read_stars('snow.grb')
snow.stars %>% as_tibble() -> snow.stars.tibble
snow.stars.tibble %>% summary(mean=mean(snow.grb, na.rm=T), min=min(snow.grb, na.rm=T), max=max(snow.grb, na.rm=T), sd=sd(snow.grb, na.rm=T))    

# snow.raster = raster('snow.grb')
# cellStats(snow.raster, stat='min', na.rm=T)
# cellStats(snow.raster, stat='max', na.rm=T)

with gdalinfo the stats are (gdalinfo -mm -stats snow.grb) Minimum=-0.000, Maximum=0.910, Mean=0.043, StdDev=0.115

The min/max values returned by gdalinfo are matching with qgis as well.

with stars, I'm getting these. max: 0.917, min: -1.17e-6, sd: 0.106 and mean: 0.0384

What am I missing?

Answer 1

You may have an XML metadata file created from a different version of your file.

If I start with just the grib file:

$ ls -l 
total 192
-rw-r--r-- 1 rowlings rowlings 189591 Apr 12 11:53 snow.grb
$

Then the stats are okay:

$  gdalinfo -mm -stats snow.grb 

Min=-0.000 Max=0.917   Computed Min/Max=-0.000,0.917
Minimum=-0.000, Maximum=0.917, Mean=0.038, StdDev=0.106
    STATISTICS_MAXIMUM=0.91749882698059
    STATISTICS_MEAN=0.038367013261481
    STATISTICS_MINIMUM=-1.1740423815354e-06
    STATISTICS_STDDEV=0.10553417638978

But GDAL has created an XML file where it keeps them:

$ ls
snow.grb  snow.grb.aux.xml

If I edit that file and set the mean to 196, then running gdalinfo again, even with the -mm option, it gets the 196 from the XML:

$ gdalinfo -stats -mm snow.grb 
  Min=-0.000 Max=0.917   Computed Min/Max=-0.000,0.917
  Minimum=-0.000, Maximum=0.917, Mean=196.000, StdDev=0.106
    STATISTICS_MAXIMUM=0.91749882698059
    STATISTICS_MEAN=196
    STATISTICS_MINIMUM=-1.1740423815354e-06
    STATISTICS_STDDEV=0.10553417638978

I can't see any option to gdalinfo to force recomputation of the metadata XML file, so you probably have to delete it.