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');


Monday, December 4, 2017

DVD, Core Meteorology

A series of DVD, director: Ron Meyer. Runtime: 30 minutes.
For grade 7 - College


It is the atmosphere that creates relatively constant daytime versus nighttime temperatures, like the temperatures found on the earth. In our solar system, no atmosphere means no possibility of life.
We live on the earth. More precisely, we live in the troposphere of the earth, like fish in the water.
The stratosphere is home to the jet stream. Also home to Ozone.
The atmosphere is a delicate dynamic balance. A balance absolutely critical for all life on the planet.
heat conduction: the pan doesn’t move.
Remarkably, the most important part of this story for us is that the troposphere is heated from the bottom to the top. In spite of being the furthest from the sun, the air closest to the ground is warmer than the air 5 miles closer to the sun.
The heat move around the atmosphere is called weather. All weather happens in the lower atmosphere.


3 factor affect weather prediction:
  1. initial state has some gap information
  2. model is not perfect
  3. the principle of time and size


Climate is a pattern over time.
Human is adapted to these stable climate patterns as well.

NG- Natural disasters

1993 Storm of the Century

3 different weather models predict different storm track.
No matter how well we can predict this, we can’t stop it from occurring. If we had control of it, the meteorologist would make more money than professional sports players make.
I may lose a few toes, not a big deal… you are in good spirit.

2004 Indian Ocean earthquake and tsunami

Earthquake is actually good to release earth’s internal energy. A solution is to create a small man-made earthquake to defuse a big one.
Experimenting with mother nature is a very difficult thing. We may trigger a bigger earthquake than the one we were trying to prevent.
Inducing quakes is an inexact science.
A 30 billion question: how often does a city like New Orleans get hit by a Category Five hurricane? Experiments show New Orleans is living in its borrowing time. A major hurricane is overdue.
Project Stormfury is trying to intervene the hurricane, but failed.
The benefit of a hurricane is redistributing heat in the atmosphere.

Tornado intercept

Tim Samra died in 2013 in Oklahoma, the first known death as a storm chaser.
Tim Samaras did not seem startled by the question from his love child, Matt Winter. “Matt,” he replied, “Kathy’s a strong woman. She understands this is my passion. And if something happened to me, she’d move on.”