2

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?

  • Can you edit your question to show the numbers that you see? You seem to have left them out. I don't see any difference. – Spacedman Apr 12 at 7:32
  • Edited the question with the numbers. – Surya Apr 12 at 10:51
  • I tried with gdal 2.3.3 and 3.0.2. Both give same result. with .grb and .tif (converted using gdal_translate) – Surya Apr 12 at 11:22
  • I think i got the problem. The gdalinfo was using the .aux.xml file. It must have contained the approximate stats. I removed the xml file and gdalinfo results are matching with stars. Thanks much for your help – Surya Apr 12 at 11:29
3

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.

| improve this answer | |
  • Thanks. I deleted the xml file and the results are matching. – Surya Apr 12 at 11:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.