Sunday, December 24, 2017

Read grb2 file in Shell, Python,R and Matlab


GRIB, or GRIdded Binary, is a concise data format by the World Meteorological Organization to store weather data.
GRIB files are a collection of self-contained records of 2D data, and the individual records stand alone as meaningful data, with no references to other records or to an overall schema. So collections of GRIB records can be appended to each other or the records separated.
The latest standard is GRIB2.
Why GRIB2?
  • compression rate. File size is 2 or 3 times smaller than NetCDF file.
  • Models are getting bigger faster than the bandwidth increases.
  • Packing is transparent to the user: self-explanatory.
Problems with GRIB: No way in GRIB to describe a collection of GRIB records.
  • Each record is independent, with no way to reference the GRIB writer’s intended schema
  • No foolproof way to combine records into the multidimensional arrays from which they were derived.
Because of such problems, any Python/R/Matlab API or Panoply reading grb2 file will lose the metadata of the records and each record is comined in a random way. Consequently, multiple records are blended with each other and become useless. So the most reliable way is still use wgrb2 in shell or embedbed shell command in C++ to extract the exact records you need in your intended order.

shell: wgrib2 C library

port install will automatically install 24 dependecies including glib2,hdf5, jasper, jpeg, netcdf, openssl, proj, zlib.
Download macport from and install.
sudo port install wgrib2
wgrib2 gfs_4_20170830_1800_024.grb2  # print out all inv rows
wgrib2 gfs_4_20170830_1800_024.grb2 -s | grep ":PRES:surface" | wgrib2 -i gfs_4_20170830_1800_024.grb2 -text data1.txt  # output surface pressure
wgrib2 gfs_4_20170830_1800_024.grb2 -s | grep ":PRES:surface" | wgrib2 -i gfs_4_20170830_1800_024.grb2 -netcdf  #
wgrib trick:


This module is a python interface to the GRIB API C library from the European Centre for Medium-Range Weather Forecasts (ECMWF)
conda install -c conda-forge pygrib
import pygrib
# Library not loaded: @rpath/libpng16.16.dylib
# Reason: Incompatible library version: requires version 51.0.0 or later, but libpng16.16.dylib provides version 49.0.0


python 3 is not compatible with pynio.
conda create --name=py2 python=2.7 anaconda
conda install -c conda-forge pynio
still incompatible
import Nio
# dlopen(/Users/jiang/anaconda3/envs/py2/lib/python2.7/site-packages/PyNIO/, 2): Library not loaded: @rpath/libgssapi_krb5.2.2.dylib
# Reason: image not found


The ECMWF GRIB-API is an application program interface accessible from C, FORTRAN and Python programs developed for encoding
Please note that GRIB-API support is being discontinued at the end of 2018.


I couldn’t find a python library to call it.

R: gdal

The Geospatial Data Abstraction Library (GDAL) is a computer software library) for reading and writing raster and vector geospatial data formats,
folder = "/Users/jiang/data/goesr/20170831CONUS"
file = "gfs_4_20170830_1800_024.grb2"
grib <- readGDAL(file)  # Large SaptialGridDataFrame, 827 Mb
# however, this is not useful because the extracted records are not in the same order with the .inv file
extract grb2 file to nc file, unpack and plot
folder = "/Users/jiang/data/gfs"
file = "gfs_4_20170830_1800_024.grb2"
# system("wgrib2 -s gfs_4_20170830_1800_024.grb2 | grep :TMP:  | wgrib2 -i gfs_4_20170830_1800_024.grb2 -netcdf",intern=T) # run once or run in shell
nc_version()  # "ncdf4_1.16_20170401"
temp <- nc_open("")
t2m.mean <- ncvar_get(temp,"TMP_2maboveground")
x <- ncvar_get(temp,"longitude")  #720
y <- ncvar_get(temp,"latitude")   #361
# dim(temp)  # 720,361
# c(nrow(temp),ncol(temp)) # 720, 361
rgb.palette <- colorRampPalette(c("snow1","snow2","snow3","seagreen","orange","firebrick"), space = "rgb")#colors
image.plot(x,y,t2m.mean,col=rgb.palette(200),main=as.expression(paste("GFS 24hr Average 2M Temperature",day,"00 UTC",sep="")),axes=T,legend.lab="o C")
plot(wrld_simpl, add = TRUE)


It actually uses 2 free toolboxes:
  1. NCtoolbox
  2. M_MAP
url = ['', mydate,'/nww3',mydate,'_00z'];
nco=ncgeodataset(url);  % Instantiate the data set
sort(nco.variables);  % list variables
waveheight=nco{'htsgwsfc'}(1,:,:);% surface 
% significant height of wind waves
waveheight = squeeze(waveheight); % squeeze shape

% Plot the field using M_MAP.  Start with setting the map
% projection using the limits of the lat/lon data itself:
m_proj('miller','lat',[min(lat(:)) max(lat(:))],'lon',[min(lon(:)) max(lon(:))]);
% Next, plot the field using the M_MAP version of pcolor.
shading flat;
% Add a coastline and axis values.
m_coast('patch',[.7 .7 .7]);
% Add a colorbar and title.
title('Example 1: WAVEWATCH III Significant Wave Height from NOMADS');