FpInputMetGribinfo: grib_domain.py

File grib_domain.py, 2.5 KB (added by pesei, 3 years ago)
Line 
1#!/usr/bin/python
2# -*- coding: utf-8 -*-
3"""
4Created on Mar 7, 2016
5
6@author: petra
7
8Version 2: return also vertical levels, based on U wind
9
10Usage: `./grib_domain.py <gribfile> [map]`
11
12The map will only be shown if the `map` keyward is present.
13
14"""
15
16import pygrib 
17import sys
18
19from mpl_toolkits.basemap import Basemap
20import numpy as np
21import matplotlib.pyplot as plt
22#from matplotlib.backends.backend_pdf import PdfPages
23 
24if len(sys.argv)==1:
25  fname='/home/petra/P/Flexpart/FpTests/Hasod_ec_hires/ER12101512'
26elif len(sys.argv)==2:
27 if sys.argv[1] == '-h' or sys.argv[1] == '--help':
28   print '\nUsage: ./grib_domain.py <gribfile> [map]'
29   print 'The map will only be shown if the `map` keyword is present.\n'
30   quit()
31else:
32  fname=sys.argv[1]
33 
34Map=False
35if len(sys.argv)==3:
36  if sys.argv[2]=='map': Map=True 
37
38en=pygrib.open(fname)
39lats, lons = en[1].latlons()
40ymin=lats.min()
41ymax=lats.max()
42xmin=lons.min()
43xmax=lons.max()
44
45ym = 0.5*(ymin+ymax)
46xm = 0.5*(xmin+xmax)
47numx=lats.shape[1]-1
48numy=lats.shape[0]-1
49dx=(xmax-xmin)/numx
50dy=(ymax-ymin)/numy
51levels=[]
52for field in en:
53  content = field.__repr__()
54  print field
55  if 'U component of wind' in content:
56    levels.append(content.split(':')[5].split(' ')[1])
57if len(levels) == 0:
58  print fname+':','x=',xmin,xmax,dx,numx+1, '  y=',ymin,ymax,dy,numy+1,\
59  'nuvz=',len(levels)
60else:
61  print fname+':','x=',xmin,xmax,dx,numx+1, '  y=',ymin,ymax,dy,numy+1,\
62  'nuvz=',len(levels),'(',min(levels),max(levels),')'
63
64if Map==True:
65
66  # lon_0 is central longitude of projection.
67  # resolution = 'c' means use crude resolution coastlines.
68  #m = Basemap(projection='robin',lon_0=45.,resolution='c')
69  #m = Basemap(projection='ortho',lat_0=ym,lon_0=xm,resolution='l')
70  xmin5=max(xmin-5.,-180.)
71  xmax5=min(xmax+5.,180.)
72  ymin5=max(ymin-5.,-90.)
73  ymax5=min(ymax+5.,90.) 
74  xmin10=max(xmin-10.,-180.)
75  xmax10=min(xmax+10.,180.)
76  ymin10=max(ymin-10.,-90.)
77  ymax10=min(ymax+10.,90.) 
78
79  m = Basemap(llcrnrlon=xmin5,llcrnrlat=ymin5,
80              urcrnrlon=xmax5,urcrnrlat=ymax5,
81              projection='mill',lat_0=ym,lon_0=xm,
82              resolution ='l',area_thresh=10000.)
83  m.drawcoastlines(linewidth=.5)
84  m.drawcountries(color='w')
85  m.fillcontinents()
86  # draw parallels and meridians.
87  m.drawparallels(np.arange(ymin10,ymax10,5.),linewidth=1.,labels=[1,0,0,0])
88  m.drawmeridians(np.arange(xmin10,xmax10,5.),linewidth=1.,latmax=90.,
89                  labels=[0,0,0,1])
90  m.plot(lons,lats,'b+',ms=2.,latlon=True)
91  plt.title('DX='+str(dx)+'  DY='+str(dy),color='b')
92  #m.etopo(scale=0.25)
93  plt.show()
94  print 'done'
hosted by ZAMG