movie_nemo_globe.py 6.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234
  1. #!/usr/bin/env python
  2. # B a r a K u d a
  3. #
  4. #
  5. # L. Brodeau, 2017
  6. import sys
  7. import os
  8. from string import replace
  9. import numpy as nmp
  10. from netCDF4 import Dataset
  11. from mpl_toolkits.basemap import Basemap
  12. #from mpl_toolkits.basemap import shiftgrid
  13. import matplotlib as mpl
  14. mpl.use('Agg')
  15. import matplotlib.pyplot as plt
  16. import matplotlib.colors as colors
  17. import warnings
  18. warnings.filterwarnings("ignore")
  19. import datetime
  20. import barakuda_colmap as bcm
  21. import barakuda_tool as bt
  22. #nsts = 4 ; # sub time-steping for creating interpolated frames (smoother video)
  23. nsts = 1 ; # sub time-steping for creating interpolated frames (smoother video)
  24. long_start = 0 ; # longitude to start the movie from...
  25. latitude = 15. ; # fixed latitude to view from
  26. lforce_mask = True
  27. year_ref_ini = 1990
  28. jt0 = 0
  29. jtN = 0
  30. #jt0 = 0
  31. #CTATM = 'T255'
  32. CTATM = 'T1279'
  33. cbox = 'Pac'
  34. fig_type='png'
  35. narg = len(sys.argv)
  36. if narg < 6: print 'Usage: '+sys.argv[0]+' <file> <variable> <LSM_file> <first record> <last record>'; sys.exit(0)
  37. cf_in = sys.argv[1] ; cv_in=sys.argv[2] ; cf_lsm=sys.argv[3]
  38. jt0 = int(sys.argv[4])
  39. jtN = int(sys.argv[5])
  40. roffset = 0.
  41. if cv_in == 'curl_ssu':
  42. cfield = 'RV(SSU)'
  43. tmin=-40. ; tmax=40. ; dtemp = 5.
  44. cpal_fld = 'ncview_blue_red'
  45. cunit = r'$[10^{-6}s^{-1}]$'
  46. cb_jump = 1
  47. if cv_in == 'mod_ssu':
  48. cfield = 'Surf. Current'
  49. tmin=0. ; tmax=3. ; dtemp = 0.25
  50. cpal_fld = 'ncview_hotres'
  51. cunit = r'$[m/s]$'
  52. cb_jump = 1
  53. if cv_in == 'sosstsst':
  54. cfield = 'SST'
  55. tmin=-2. ; tmax=32. ; dtemp = 1.
  56. cpal_fld = 'ncview_nrl'
  57. cunit = r'$^{\circ}C$'
  58. cb_jump = 2
  59. elif cv_in == 'somxl010':
  60. cfield == 'MLD'
  61. tmin=50. ; tmax=1500. ; dtemp = 50.
  62. cpal_fld = 'viridis_r'
  63. #bt.chck4f(cf_ice)
  64. bt.chck4f(cf_in)
  65. id_fld = Dataset(cf_in)
  66. vtime = id_fld.variables['time_counter'][:]
  67. id_fld.close()
  68. Nt = len(vtime)
  69. if jtN > Nt: print 'ERROR: your jtN > Nt !!!'; sys.exit(0)
  70. if jtN == 0: jtN = Nt
  71. bt.chck4f(cf_lsm)
  72. id_lsm = Dataset(cf_lsm)
  73. #LOLO for curl => F !!!
  74. if lforce_mask: XMSK = id_lsm.variables['fmask'][0,0,:,:]
  75. XLON = id_lsm.variables['glamf'][0,:,:]
  76. XLAT = id_lsm.variables['gphif'][0,:,:]
  77. id_lsm.close()
  78. (nj,ni) = nmp.shape(XLON)
  79. #pmsk = nmp.ma.masked_where(XMSK[:,:] > 0.2, XMSK[:,:]*0.+40.)
  80. #idx_oce = nmp.where(XMSK[:,:] > 0.5)
  81. XFLD = nmp.zeros((nj,ni))
  82. if nsts>1:
  83. XFLDt = nmp.zeros((nj,ni))
  84. XFLDtp1 = nmp.zeros((nj,ni))
  85. params = { 'font.family':'Ubuntu',
  86. 'font.size': int(12),
  87. 'legend.fontsize': int(12),
  88. 'xtick.labelsize': int(12),
  89. 'ytick.labelsize': int(12),
  90. 'axes.labelsize': int(12) }
  91. mpl.rcParams.update(params)
  92. cfont_clb = { 'fontname':'Arial', 'fontweight':'normal', 'fontsize':13, 'color':'white' }
  93. cfont_title = { 'fontname':'Ubuntu Mono', 'fontweight':'normal', 'fontsize':18, 'color':'white' }
  94. cfont_mail = { 'fontname':'Times New Roman', 'fontweight':'normal', 'fontstyle':'italic', 'fontsize':10, 'color':'0.5'}
  95. # Colormaps for fields:
  96. pal_fld = bcm.chose_colmap(cpal_fld)
  97. norm_fld = colors.Normalize(vmin = tmin, vmax = tmax, clip = False)
  98. #pal_ice = bcm.chose_colmap(cpal_ice)
  99. #norm_ice = colors.Normalize(vmin = rmin_ice, vmax = 1, clip = False)
  100. pal_lsm = bcm.chose_colmap('blk')
  101. norm_lsm = colors.Normalize(vmin = 0., vmax = 1., clip = False)
  102. vc_fld = nmp.arange(tmin, tmax + dtemp, dtemp)
  103. print ''
  104. rh = 7.5
  105. jrot = -1 + jt0*nsts
  106. print '\n jt0, jtN =>', jt0, jtN, '\n'
  107. for jt in range(jt0,jtN):
  108. ct = '%3.3i'%(jt+1)
  109. cd = str(datetime.datetime.strptime('1990 '+ct, '%Y %j'))
  110. cdate = cd[:10] ; print ' *** Date :', cdate
  111. print "\n *** Reading record #"+str(ct)+" of "+cv_in+" in "+cf_in
  112. id_fld = Dataset(cf_in)
  113. if nsts==1:
  114. XFLD = id_fld.variables[cv_in][jt,:,:]
  115. else:
  116. if jt > jt0:
  117. XFLDt[:,:] = XFLDtp1[:,:]
  118. else:
  119. XFLDt = id_fld.variables[cv_in][jt,:,:]
  120. XFLDtp1 = id_fld.variables[cv_in][jt+1,:,:]
  121. id_fld.close()
  122. print " => done!"
  123. # sub time - step
  124. if nsts>1: dFdt = (XFLDtp1[:,:] - XFLDt[:,:])/float(nsts)
  125. for js in range(nsts):
  126. jrot = jrot+1
  127. # Linear interpolation: # yN = y1 + (tN-t1)*(y2-y1)/(t2-t1)
  128. if nsts>1: XFLD[:,:] = XFLDt[:,:] + float(js)*dFdt
  129. rot = (long_start + (0.5/float(nsts)*float(jrot)))%360.
  130. rot = -rot
  131. cfig = 'figs/'+cv_in+'_NEMO'+'_d'+ct+'_'+str(js)+'.'+fig_type
  132. print ' *** reference longitude =', rot
  133. fig = plt.figure(num = 1, figsize=(rh,1.167*rh), dpi=None, facecolor='b', edgecolor='k')
  134. ax = plt.axes([0.005, 0.05, 0.99, 0.99], facecolor = '0.35')
  135. plt.title('Ocean (NEMO@ORCA12 + IFS@'+CTATM+'): '+cfield+', '+cdate, **cfont_title)
  136. print ' *** Creating new projection'
  137. carte = Basemap(projection='ortho', lat_0=latitude, lon_0=rot, resolution='h')
  138. x0,y0 = carte(XLON,XLAT)
  139. if lforce_mask: XFLD = nmp.ma.masked_where(XMSK[:,:] < 0.5, XFLD[:,:])
  140. print ' *** Ploting on map...'
  141. cf = carte.pcolor(x0, y0, XFLD, cmap=pal_fld, norm=norm_fld)
  142. if lforce_mask: cc = carte.contour(x0, y0, XMSK, [ 0. ], colors='k', linewidths=0.75)
  143. plt.annotate('L. Brodeau, brodeau@gmail.com', xy=(0.7, 0.1), xycoords='figure fraction', **cfont_mail)
  144. # Colorbar:
  145. ax2 = plt.axes([0.005, 0.06, 0.99, 0.025])
  146. clb = mpl.colorbar.ColorbarBase(ax2, ticks=vc_fld, cmap=pal_fld, norm=norm_fld, orientation='horizontal', extend='both')
  147. if cb_jump > 1:
  148. cb_labs = [] ; cpt = 0
  149. for rr in vc_fld:
  150. if cpt % cb_jump == 0:
  151. cb_labs.append(str(int(rr)))
  152. else:
  153. cb_labs.append(' ')
  154. cpt = cpt + 1
  155. clb.ax.set_xticklabels(cb_labs)
  156. clb.set_label(cunit, **cfont_clb)
  157. clb.ax.yaxis.set_tick_params(color='white') ; # set colorbar tick color
  158. clb.outline.set_edgecolor('white') ; # set colorbar edgecolor
  159. plt.setp(plt.getp(clb.ax.axes, 'xticklabels'), color='white') ; # set colorbar ticklabels
  160. print ' *** Saving figure...'
  161. #plt.savefig(cfig, dpi=160, orientation='portrait', facecolor='#06051F')
  162. plt.savefig(cfig, dpi=160, orientation='portrait', facecolor='k')
  163. print ' => '+cfig+' created!\n'
  164. plt.close(1)
  165. del cf, carte, x0, y0
  166. if nsts>1: del dFdt