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

Last change on this file since 5066 was 4700, checked in by musat, 14 months ago

Adaptation scripts AXE4 a spirit2
IonelaMusat?

  • Property svn:executable set to *
File size: 4.4 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  inputfile = ''
20  outputdir = ''
21  errmsg='Use: '+str(sys.argv[0])+' -i <inputfile> -o <outputdir>'
22  try:
23     opts, args = getopt.getopt(argv,"h:i:o:")
24  except getopt.GetoptError:
25     print(errmsg)
26     sys.exit(2)
27  for opt, arg in opts:
28     if opt == '-h':
29        print(errmsg)
30        sys.exit()
31     elif opt in ("-i"):
32        inputfile = arg
33     elif opt in ("-o"):
34        outputdir = arg
35  if len(inputfile) == 0 or len(outputdir) == 0:
36    print('Please specify an input file and an output directory.')
37    print(errmsg)
38    sys.exit()
39  print('Input file is ', inputfile)
40  print('Output directory is ', outputdir)
41
42  # MAIN PROGRAM
43  # -----------------------------------------------------------------
44 
45  ncfile=inputfile
46  shortname=(ncfile.rsplit('/',1))[1].rsplit('_',-1)[0].replace(".","-")
47 
48  # LOADING THE GCM RESULTS
49  # -----------------------------------------------------------------
50 
51  ncgcm = netCDF4.Dataset(ncfile)
52  longi0 = ncgcm.variables['lon'][:] 
53  longi = longi0
54  lati0 = ncgcm.variables['lat'][:]
55  lati = lati0
56  snowgcm = ncgcm.variables['snow'][:] # precip rate in mm/hr
57 
58  ncpath='/thredds/ipsl/fabric/lmdz/AXE4/'
59  ncfractlic=ncpath + 'fract_lic_144x142.nc'
60  ncmask = netCDF4.Dataset(ncfractlic)
61  fract_lic = ncmask.variables['fract_lic'][:]
62  fract_lic = fract_lic[0,:,:] # constant anyway
63 
64  # DISPLAYING THE RESULTS
65  # -----------------------------------------------------------------
66 
67  # POLAR MAP
68  # source: http://polar.ncep.noaa.gov/global/examples/usingpython.shtml
69  #         http://matplotlib.org/basemap/users/pstere.html
70  # -----------------------------------------------------------------
71 
72  # Some postprocessing
73  #>>> snow.shape
74  #(1, 143, 144)
75  # Convert missing values to nan
76  snowgcm = np.nanmean(snowgcm,axis=0)*86400.*365.
77  snowgcm[fract_lic==0]=np.nan
78  data = snowgcm
79 
80  # Setup north polar stereographic basemap.
81  # The longitude lon_0 is at 6-o'clock, and the
82  # latitude circle boundinglat is tangent to the edge
83  # of the map at lon_0. Default value of lat_ts
84  # (latitude of true scale) is pole.
85  m = Basemap(projection='spstere',boundinglat=-60, \
86    lon_0=180.,resolution='c')
87 
88  # Convert the lat/lon values to x/y projections.
89 
90  Lon, Lat = meshgrid(longi,lati)
91  x, y = m(Lon,Lat)
92 
93  # draw filled contours.
94  clevs = [0.,20,50,100,200,300,400,500,700,1E3,1500,2E3,7E3]
95  #clevs = [0.,0.25,0.5,0.75,1.,1.25,1.5,2.,3.,5.,10.,20]
96  cs = m.contourf(x,y,data,clevs, \
97  # norm = mpl.colors.PowerNorm(gamma=0.15), \
98  # norm = mpl.colors.LogNorm(vmin=clevs[0],vmax=clevs[len(clevs)-1]), \
99  # norm = mpl.colors.Normalize(vmin=clevs[0],vmax=clevs[len(clevs)-1]), \
100   norm = mpl.colors.BoundaryNorm(boundaries=clevs, ncolors=256), \
101   cmap='jet')
102  #cs = m.pcolor(x,y,data, \
103  #  norm = mpl.colors.Normalize(vmin=clevs[0],vmax=clevs[len(clevs)-1]))
104  #cs = m.contourf(x,y,data,clevs)
105  #cs = m.contour(x,y,data,clevs,colors='k',linewidths=1.)
106  # Add a colorbar and title, and number of bins
107  cbar = m.colorbar(cs,location='bottom',pad="5%")
108  #cbar.set_label('mm')
109  #cbar.locator = ticker.MaxNLocator(nbins=15)
110  #cbar.update_ticks()
111 
112  # Next, plot the field using the fast pcolormesh routine and set the colormap to jet.
113  #cs = m.pcolormesh(x,y,data,shading='flat', \
114  #cmap=plt.cm.jet, \
115  #vmin=0., vmax=100)
116 
117  # Add a coastline and axis values.
118  m.drawcoastlines()
119  #m.fillcontinents()
120  m.drawmapboundary()
121  m.drawparallels(np.arange(-90.,95.,5.), \
122  labels=[1,0,0,0])
123  m.drawmeridians(np.arange(-180.,180.,30.), \
124  labels=[0,0,0,1])
125 
126  plt.title('Precip in mm/yr ('+shortname+')')
127 
128  savefig(outputdir+'/'+'precipSH-'+shortname+'.png', bbox_inches='tight')
129  #plt.show()
130 
131  # -----------------------------------------------------------------
132  # -----------------------------------------------------------------
133  # -----------------------------------------------------------------
134
135# MAIN PROGRAM
136# Just calling the main function
137if __name__ == "__main__":
138   main(sys.argv[1:])
Note: See TracBrowser for help on using the repository browser.