source: trunk/src/scripts_Laura/XCPGR_SSMIS_test.py @ 45

#!/usr/bin/env python
# -*- coding: utf-8 -*-
import string
import numpy as np
import matplotlib.pyplot as plt
from pylab import *
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.basemap import shiftgrid, cm
import netCDF4
import ffgrid2





########################
# EVOLUTION TEMPORELLE #
########################
## Autour d une zone de glace de mer (lon entre -70 et -90 ; lat entre -70 et -80) - mois: JUNE ##
## ch12 ##
bbzone_CH12_JUN = nonzero((lon12_JUN >= -90.) & (lon12_JUN <= -70.) & (lat12_JUN >= -80.) & (lat12_JUN <= -70.))
tb12_JUN_moy = np.zeros([30],float)

## ch13 ##
bbzone_CH13_JUN = nonzero((lon13_JUN >= -90.) & (lon13_JUN <= -70.) & (lat13_JUN >= -80.) & (lat13_JUN <= -70.))
tb13_JUN_moy = np.zeros([30],float)

## ch15 ##
bbzone_CH15_JUN = nonzero((lon15_JUN >= -90.) & (lon15_JUN <= -70.) & (lat15_JUN >= -80.) & (lat15_JUN <= -70.))
tb15_JUN_moy = np.zeros([30],float)

## ch 16 ##
bbzone_CH16_JUN = nonzero((lon16_JUN >= -90.) & (lon16_JUN <= -70.) & (lat16_JUN >= -80.) & (lat16_JUN <= -70.))
tb16_JUN_moy = np.zeros([30],float)

for ijr in range (0,30):
     print 'jour=', ijr+1
     ## ch12 ##
     ind_jr12 = np.where(jjr12_JUN[bbzone_CH12_JUN]==ijr+1)[0]
     a = tb12_JUN[bbzone_CH12_JUN][ind_jr12]
     tb12_JUN_moy[ijr] = mean(a[nonzero((a!=-500.)&(a!=0.))])
     ## ch13 ##
     ind_jr13 = np.where(jjr13_JUN[bbzone_CH13_JUN]==ijr+1)[0]
     b = tb13_JUN[bbzone_CH13_JUN][ind_jr13]
     tb13_JUN_moy[ijr] = mean(b[nonzero((b != -500.) & (b != 0.))])
     ## ch15 ##
     ind_jr15 = np.where(jjr15_JUN[bbzone_CH15_JUN]==ijr+1)[0]
     c = tb15_JUN[bbzone_CH15_JUN][ind_jr15]
     tb15_JUN_moy[ijr] = mean(c[nonzero((c != -500.) & (c != 0.))])
     ## ch16 ##
     ind_jr16 = np.where(jjr16_JUN[bbzone_CH16_JUN]==ijr+1)[0]
     d = tb16_JUN[bbzone_CH16_JUN][ind_jr16]
     tb16_JUN_moy[ijr] = mean(d[nonzero((d != -500.) & (d != 0.))])


XCPGR = (tb12_JUN_moy - tb16_JUN_moy) / (tb12_JUN_moy + tb16_JUN_moy)

## ch12 - ch13 // ch15 - ch16  ##
fig=plt.figure()
plt.subplot(2,1,1)
plt.plot(arange(0,30,1), tb12_JUN_moy, label = '19.35GHz-H', c = 'r')
plt.plot(arange(0,30,1), tb13_JUN_moy, label = '19.35GHz-V', c = 'b')
plt.plot(arange(0,30,1), tb15_JUN_moy, label = '37GHz-H', c = 'm')
plt.plot(arange(0,30,1), tb16_JUN_moy, label = '37GHz-V', c = 'g')
plt.xticks(arange(0,30,1), arange(1,31,1))
grid(True)
ylabel('brightness temperature')
legend(loc=7)
plt.title('SSMIS - JUNE2010')
plt.subplot(2,1,2)
plt.plot(arange(0,30,1), XCPGR)
plt.xticks(arange(0,30,1), arange(1,31,1))
grid(True)
ylabel('Cross polarized gradient ratio')
fig.show()


################
# CARTOGRAPHIE #
################
## Etude sur l'Antarctique ##
dx=5.
dy=5.
x0, x1 = -180, 180
y0, y1 = -90, -30

## ch12 ##
bbtb_ch12_JUN=nonzero((tb12_JUN!=-500.)&(tb12_JUN!=0.))
OUTZCH12_JUN = np.zeros([len(np.arange(y0, y1+1, dy)),len(np.arange(x0, x1+1, dx)),30], float)
outzch12_JUN = np.zeros([len(np.arange(y0, y1+1, dy)),len(np.arange(x0, x1+1, dx))], float)
lonch12_JUN=np.zeros([len(np.arange(x0, x1+1, dx))], float)
latch12_JUN=np.zeros([len(np.arange(y0, y1+1, dy))], float)

## ch16 ##
bbtb_ch16_JUN=nonzero((tb16_JUN!=-500.)&(tb16_JUN!=0.))
OUTZCH16_JUN = np.zeros([len(np.arange(y0, y1+1, dy)),len(np.arange(x0, x1+1, dx)),30], float)
outzch16_JUN = np.zeros([len(np.arange(y0, y1+1, dy)),len(np.arange(x0, x1+1, dx))], float)
lonch16_JUN=np.zeros([len(np.arange(x0, x1+1, dx))], float)
latch16_JUN=np.zeros([len(np.arange(y0, y1+1, dy))], float)

for ijr in range(0,30):
    print 'jour=', ijr+1
    ## ch12 ##
    ind_jr12_JUN = np.where(jjr12_JUN[bbtb_ch12_JUN] == ijr+1)[0]
    xx = lon12_JUN[bbtb_ch12_JUN][ind_jr12_JUN]
    yy = lat12_JUN[bbtb_ch12_JUN][ind_jr12_JUN]
    zz = tb12_JUN[bbtb_ch12_JUN][ind_jr12_JUN]
    zz0 = min(zz)
    zz1 = max(zz)
    outz, outx, outy = ffgrid2.ffgrid(xx, yy, zz, dx, dy, x0,x1,y0,y1,zz0, zz1)
    outzch12_JUN=outz
    lonch12_JUN=outx
    latch12_JUN=outy
    OUTZCH12_JUN[:,:,ijr] = outzch12_JUN[:,:]
    ## ch16 ##
    ind_jr16_JUN = np.where(jjr16_JUN[bbtb_ch16_JUN] == ijr+1)[0]
    xx = lon16_JUN[bbtb_ch16_JUN][ind_jr16_JUN]
    yy = lat16_JUN[bbtb_ch16_JUN][ind_jr16_JUN]
    zz = tb16_JUN[bbtb_ch16_JUN][ind_jr16_JUN]
    zz0 = min(zz)
    zz1 = max(zz)
    outz, outx, outy = ffgrid2.ffgrid(xx, yy, zz, dx, dy, x0,x1,y0,y1,zz0, zz1)
    outzch16_JUN=outz
    lonch16_JUN=outx
    latch16_JUN=outy
    OUTZCH16_JUN[:,:,ijr] = outzch16_JUN[:,:]


XCPGR = (OUTZCH12_JUN - OUTZCH16_JUN) / (OUTZCH12_JUN + OUTZCH16_JUN)

for ijr in range (13,20):
    figure()     
    plt.ion()
    m = Basemap(llcrnrlon=-180, urcrnrlon=180, llcrnrlat=-90, urcrnrlat=-30, projection='cyl', resolution='c', fix_aspect=True)
    m.drawcoastlines(linewidth = 1)
    m.drawparallels(np.arange(-90., 90., 20))
    m.drawmeridians(np.arange(-180., 180., 20))
    #m.fillcontinents()
    xii,yii = m(*np.meshgrid(lonch12_JUN, latch12_JUN))
    clevs = arange(-0.3, 0.01, 0.01)
    cs = m.contourf(xii, yii, XCPGR[:,:,ijr], clevs, cmap=cm.s3pcpn_l_r)
    cbar = colorbar(cs)
    cbar.set_label('XCPGR SSMIS - JUNE')
    

##########################
# DIAGRAMME DE HOVMOLLER #
##########################
## Etude sur une tranche de latitude (proche de la cote: entre lat0 = -80 et lat1 = -75 - mois: JUNE ##
# shape(tbch17_anom_JUN) = [ilat, ilon, ijr]
bbtranche12_JUN = nonzero((latch12_JUN >= -80.) & (latch12_JUN <= -75.))
grad_ratio = np.zeros([len(lonch12_JUN), 30], float)
for ilon in range (0, len(lonch12_JUN)):
    for ijr in range (0,30):
        grad_ratio[ilon, ijr] = XCPGR[bbtranche12_JUN][:, ilon, ijr].mean()


y_time, x_space = m(*np.meshgrid(arange(0,30,1), lonch12_JUN))
fig = plt.figure()
plt.pcolor(x_space, y_time, grad_ratio, cmap=cm.s3pcpn_l_r, vmin = -0.3, vmax = 0.1)
plt.axis([-180., 180., 0, 29])
cb = plt.colorbar()
cb.set_label('Cross polarization gradient ratio - SSMIS CH12-CH16')
plt.xticks(arange(-180.,220.,40))
plt.yticks(arange(0, 30, 1), arange(1,31,1))
plt.xlabel('longitude')
plt.ylabel('JUNE 2010')