FpInputMetGribinfo: grib_domain3.py

#!/usr/bin/python3
# -*- coding: utf-8 -*-
"""
Created on Mar 7, 2016

@author: petra

Version 2: return also vertical levels, based on U wind
Version 3: for Python3

Usage: `./grib_domain.py <gribfile> [map]`

The map will only be shown if the `map` keyward is present.

Dependencies: you need python3-grib and python3-mpltoolkits.basemap
  (Debian packages) or equivalent packages of your OS

"""

import pygrib 
import sys

from mpl_toolkits.basemap import Basemap
import numpy as np
import matplotlib.pyplot as plt
#from matplotlib.backends.backend_pdf import PdfPages

 
if len(sys.argv)==1:
 if sys.argv[1] == '-h' or sys.argv[1] == '--help':
   print('\nUsage: ./grib_domain.py <gribfile> [map]')
   print('The map will only be shown if the `map` keyword is present.\n')
   quit()
elif len(sys.argv)>=2:
 if sys.argv[1] == '-h' or sys.argv[1] == '--help':
   print('\nUsage: ./grib_domain.py <gribfile> [map]')
   print('The map will only be shown if the `map` keyword is present.\n')
   quit()
 else:
  fname=sys.argv[1]
  
Map=False
if len(sys.argv)==3:
  if sys.argv[2]=='map': Map=True  

en=pygrib.open(fname)
lats, lons = en[1].latlons()
ymin=lats.min()
ymax=lats.max()
xmin=lons.min()
xmax=lons.max()

ym = 0.5*(ymin+ymax)
xm = 0.5*(xmin+xmax)
numx=lats.shape[1]-1
numy=lats.shape[0]-1
dx=(xmax-xmin)/numx
dy=(ymax-ymin)/numy
levels=[]
print('Fields found in',fname,':')
for field in en:
  content = field.__repr__()
  print(field)
  if 'U component of wind' in content:
    levels.append(content.split(':')[5].split(' ')[1])
if len(levels) == 0:
  print(fname+': ','x=',xmin,xmax,dx,numx+1, '  y=',ymin,ymax,dy,numy+1,\
  'nuvz=',len(levels))
else:
  print(fname+': ','x/xmin/xmax/dx/numx=',xmin,xmax,dx,numx+1, '; y=',ymin,ymax,dy,numy+1,\
  '; nuvz (min max) =',len(levels),'(',min(levels),max(levels),')')

if Map==True:

  # lon_0 is central longitude of projection.
  # resolution = 'c' means use crude resolution coastlines.
  #m = Basemap(projection='robin',lon_0=45.,resolution='c')
  #m = Basemap(projection='ortho',lat_0=ym,lon_0=xm,resolution='l')
  xmin5=max(xmin-5.,-180.)
  xmax5=min(xmax+5.,180.)
  ymin5=max(ymin-5.,-90.)
  ymax5=min(ymax+5.,90.)  
  xmin10=max(xmin-10.,-180.)
  xmax10=min(xmax+10.,180.)
  ymin10=max(ymin-10.,-90.)
  ymax10=min(ymax+10.,90.)  

  m = Basemap(llcrnrlon=xmin5,llcrnrlat=ymin5,
              urcrnrlon=xmax5,urcrnrlat=ymax5,
              projection='mill',lat_0=ym,lon_0=xm,
              resolution ='l',area_thresh=10000.)
  m.drawcoastlines(linewidth=.5)
  m.drawcountries(color='w')
  m.fillcontinents()
  # draw parallels and meridians.
  m.drawparallels(np.arange(ymin10,ymax10,5.),linewidth=1.,labels=[1,0,0,0])
  m.drawmeridians(np.arange(xmin10,xmax10,5.),linewidth=1.,latmax=90.,
                  labels=[0,0,0,1])
  m.plot(lons,lats,'b+',ms=2.,latlon=True)
  plt.title('DX='+str(dx)+'  DY='+str(dy),color='b')
  #m.etopo(scale=0.25)
  plt.show()
  print('done')