source: BOL/Multi_atlas/AXE4/make_axe4-map-precip-obs.py @ 5272

Last change on this file since 5272 was 4700, checked in by musat, 15 months ago

Adaptation scripts AXE4 a spirit2
IonelaMusat?

  • Property svn:executable set to *
File size: 4.2 KB
Line 
1#! /usr/bin/env python
2# sources:
3# http://matplotlib.org/basemap/users/examples.html
4# http://matplotlib.org/users/colormapnorms.html
5from mpl_toolkits.basemap import Basemap
6import numpy as np
7import math as m
8import matplotlib.pyplot as plt
9#import matplotlib.dates as dates
10#import datetime
11#from matplotlib.dates import date2num
12from matplotlib import ticker
13#import pandas
14from pylab import *
15import netCDF4
16import sys, getopt
17
18def main(argv):
19  outputdir = ''
20  errmsg='Use: '+str(sys.argv[0])+' -o <outputdir>'
21  try:
22     opts, args = getopt.getopt(argv,"h:o:")
23  except getopt.GetoptError:
24     print(errmsg)
25     sys.exit(2)
26  for opt, arg in opts:
27     if opt == '-h':
28        print(errmsg)
29        sys.exit()
30     elif opt in ("-o"):
31        outputdir = arg
32  if len(outputdir) == 0:
33    print('Please specify an output directory.')
34    print(errmsg)
35    sys.exit()
36  print('Output directory is ', outputdir)
37
38  # BEGINNING OF THE MAIN ROUTINE
39  # -----------------------------------------------------------------
40 
41  ncpath='/thredds/ipsl/fabric/lmdz/AXE4/'
42  ncfile=ncpath + 'CloudSat_Snowfall_rate_2007_2010_regrid.nc'
43  ncobs = netCDF4.Dataset(ncfile)
44  longi = ncobs.variables['lon'][:]
45  lati = ncobs.variables['lat'][:]
46  snowobs = ncobs.variables['Snowfall_rate'][:] # precip rate in mm/hr
47 
48  ncfractlic=ncpath + 'fract_lic_144x142.nc'
49  ncmask = netCDF4.Dataset(ncfractlic)
50  fract_lic = ncmask.variables['fract_lic'][:]
51  fract_lic = fract_lic[0,:,:] # constant anyway
52 
53  # DISPLAYING THE RESULTS
54  # -----------------------------------------------------------------
55 
56  # POLAR MAP
57  # source: http://polar.ncep.noaa.gov/global/examples/usingpython.shtml
58  #         http://matplotlib.org/basemap/users/pstere.html
59  # -----------------------------------------------------------------
60 
61  # Some postprocessing
62  #>>> snow.shape
63  #(1, 143, 144)
64  # Convert missing values to nan
65  snowobs = np.array(snowobs) # remove the empty values '--'
66  snowobs[snowobs==-9.00000000e+33]=np.nan # remove -9.00000000e+33
67  snowobs[fract_lic==0]=np.nan # only keep Antarctica or Greenland
68  data = snowobs*24.*365. # convertion to mm/yr
69 
70  # Setup north polar stereographic basemap.
71  # The longitude lon_0 is at 6-o'clock, and the
72  # latitude circle boundinglat is tangent to the edge
73  # of the map at lon_0. Default value of lat_ts
74  # (latitude of true scale) is pole.
75  m = Basemap(projection='spstere',boundinglat=-60, \
76    lon_0=180.,resolution='c')
77 
78  # Convert the lat/lon values to x/y projections.
79 
80  Lon, Lat = meshgrid(longi,lati)
81  x, y = m(Lon,Lat)
82 
83  # draw filled contours.
84  clevs = [0.,20,50,100,200,300,400,500,700,1E3,1500,2E3,7E3]
85  #clevs = [0.,0.25,0.5,0.75,1.,1.25,1.5,2.,3.,5.,10.,20]
86  cs = m.contourf(x,y,data,clevs, \
87  # norm = mpl.colors.PowerNorm(gamma=0.15), \
88  # norm = mpl.colors.LogNorm(vmin=clevs[0],vmax=clevs[len(clevs)-1]), \
89  # norm = mpl.colors.Normalize(vmin=clevs[0],vmax=clevs[len(clevs)-1]), \
90   norm = mpl.colors.BoundaryNorm(boundaries=clevs, ncolors=256), \
91   cmap='jet')
92  #cs = m.pcolor(x,y,data, \
93  #  norm = mpl.colors.Normalize(vmin=clevs[0],vmax=clevs[len(clevs)-1]))
94  #cs = m.contourf(x,y,data,clevs)
95  #cs = m.contour(x,y,data,clevs,colors='k',linewidths=1.)
96  # Add a colorbar and title, and number of bins
97  cbar = m.colorbar(cs,location='bottom',pad="5%")
98  #cbar.set_label('mm')
99  #cbar.locator = ticker.MaxNLocator(nbins=15)
100  #cbar.update_ticks()
101 
102  # Next, plot the field using the fast pcolormesh routine and set the colormap to jet.
103  #cs = m.pcolormesh(x,y,data,shading='flat', \
104  #cmap=plt.cm.jet, \
105  #vmin=0., vmax=100)
106 
107  # Add a coastline and axis values.
108  m.drawcoastlines()
109  #m.fillcontinents()
110  m.drawmapboundary()
111  m.drawparallels(np.arange(-90.,95.,5.), \
112  labels=[1,0,0,0])
113  m.drawmeridians(np.arange(-180.,180.,30.), \
114  labels=[0,0,0,1])
115 
116  plt.title('Precip in mm/yr (CloudSat obs)')
117  #plt.show()
118  savefig(outputdir+'/'+'precipSH-OBS.png', bbox_inches='tight')
119 
120  # -----------------------------------------------------------------
121  # -----------------------------------------------------------------
122  # -----------------------------------------------------------------
123
124if __name__ == "__main__":
125   main(sys.argv[1:])
126
Note: See TracBrowser for help on using the repository browser.