source: trunk/src/scripts_Laura/monthly_diff_frequ_AMSUA_test.py @ 54

#!/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
from netCDF4 import Dataset
import ffgrid2


## ATTENTION: lecture des fichiers NETCDF AMSUA - definition des variables dans "read_netcdf_AMSUA_test.py" ##
## lon
## lat
## jrs_MONTH
## chan
len_month = np.array([31, 28, 31, 30, 31, 30, 31])

###############################################
## emissivite moyenne sur la periode d etude ##
###############################################
Njr = np.sum(len_month)
emis_glob = np.zeros([len(chan), len(lat), len(lon), Njr], float)
for ich in range (0, len(chan)):
    for ilon in range (0, len(lon)):
        for ilat in range (0, len(lat)):
            e1 = np.append(emis_JAN[ich, ilat, ilon, :], emis_FEB[ich, ilat, ilon, :])
            e2 = np.append(e1, emis_MAR[ich, ilat, ilon, :])
            e3 = np.append(e2, emis_APR[ich, ilat, ilon, :])
            e4 = np.append(e3, emis_MAY[ich, ilat, ilon, :])
            e5 = np.append(e4, emis_JUN[ich, ilat, ilon, :])
            e6 = np.append(e5, emis_JUL[ich, ilat, ilon, :])
            emis_glob[ich, ilat, ilon, :] = e6


emis_glob_moy = np.zeros([len(chan), len(lat), len(lon)], float)
for ich in range (0, len(chan)):
    for ilon in range (0, len(lon)):
        for ilat in range (0, len(lat)):
            emis_glob_moy[ich, ilat, ilon] = mean(emis_glob[ich, ilat, ilon, :][nonzero(isnan(emis_glob[ich, ilat, ilon, :]) == False)])



#############################################################
## variance de l emissivite moyenne sur la periode d etude ##
#############################################################
var = np.zeros([len(chan), len(lat), len(lon)], float)
for ich in range (0, len(chan)):
    for ilon in range (0, len(lon)):
        for ilat in range (0, len(lat)):
            var[ich, ilat, ilon] = np.sum((emis_glob[ich, ilat, ilon, :][nonzero(isnan(emis_glob[ich, ilat, ilon, :]) == False)] - emis_glob_moy[ich, ilat, ilon])**2)


var_glob = var / Njr


## selection des zones en fonction de la pgade de valeur de var ##
for ich in range (0, len(chan)):
    for ilon in range (0, len(lon)):
        for ilat in range (0, len(lat)):
#            if ((var_glob[ich, ilat, ilon] >  0.0040) & (var_glob[ich, ilat, ilon] < 0.0060)):
#                continue
#            if ((var_glob[ich, ilat, ilon] >=  0.0024) & (var_glob[ich, ilat, ilon] <  0.0060)):
#                continue
            if (var_glob[ich, ilat, ilon] <=  0.0030):
                continue
            else:
                var_glob[ich, ilat, ilon] = NaN

## carto ##
plt.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., -30., 20))
m.drawmeridians(np.arange(-180., 180., 20))
xii,yii = m(*np.meshgrid(lon, lat))
clevs = (arange(0., 0.01, 0.0001))
cs = m.contourf(xii, yii, var_glob[0, :, :], clevs, cmap=cm.s3pcpn_l_r)
cbar = colorbar(cs)
cbar.set_label('CH1 - 23.8GHz')
plt.title('var(emis) - AMSUA')
xticks(arange(-180, 200, 20), rotation = 45)
yticks(arange(-90, -10, 20))



##########################################################
## biais de l emissivite moyenne sur la periode d etude ##
##########################################################

## emis moyenne par mois ##
emis_JAN_mean = np.zeros([len(chan), len(lat), len(lon)], float)
emis_FEB_mean = np.zeros([len(chan), len(lat), len(lon)], float)
emis_MAR_mean = np.zeros([len(chan), len(lat), len(lon)], float)
emis_APR_mean = np.zeros([len(chan), len(lat), len(lon)], float)
emis_MAY_mean = np.zeros([len(chan), len(lat), len(lon)], float)
emis_JUN_mean = np.zeros([len(chan), len(lat), len(lon)], float)
emis_JUL_mean = np.zeros([len(chan), len(lat), len(lon)], float)
for ich in range (0, len(chan)):
    for ilon in range (0, len(lon)):
        for ilat in range (0, len(lat)):
            emis_JAN_mean[ich, ilat, ilon] = mean(emis_JAN[ich, ilat, ilon, :][nonzero(isnan(emis_JAN[ich, ilat, ilon, :]) == False)])
            emis_FEB_mean[ich, ilat, ilon] = mean(emis_FEB[ich, ilat, ilon, :][nonzero(isnan(emis_FEB[ich, ilat, ilon, :]) == False)])
            emis_MAR_mean[ich, ilat, ilon] = mean(emis_MAR[ich, ilat, ilon, :][nonzero(isnan(emis_MAR[ich, ilat, ilon, :]) == False)])
            emis_APR_mean[ich, ilat, ilon] = mean(emis_APR[ich, ilat, ilon, :][nonzero(isnan(emis_APR[ich, ilat, ilon, :]) == False)])
            emis_MAY_mean[ich, ilat, ilon] = mean(emis_MAY[ich, ilat, ilon, :][nonzero(isnan(emis_MAY[ich, ilat, ilon, :]) == False)])
            emis_JUN_mean[ich, ilat, ilon] = mean(emis_JUN[ich, ilat, ilon, :][nonzero(isnan(emis_JUN[ich, ilat, ilon, :]) == False)])
            emis_JUL_mean[ich, ilat, ilon] = mean(emis_JUL[ich, ilat, ilon, :][nonzero(isnan(emis_JUL[ich, ilat, ilon, :]) == False)])


## biais ch1-ch15 ##
bias_emis_JAN_mean = emis_JAN_mean[0, :, :] - emis_JAN_mean[1, :, :]
bias_emis_FEB_mean = emis_FEB_mean[0, :, :] - emis_FEB_mean[1, :, :]
bias_emis_MAR_mean = emis_MAR_mean[0, :, :] - emis_MAR_mean[1, :, :]
bias_emis_APR_mean = emis_APR_mean[0, :, :] - emis_APR_mean[1, :, :]
bias_emis_MAY_mean = emis_MAY_mean[0, :, :] - emis_MAY_mean[1, :, :]
bias_emis_JUN_mean = emis_JUN_mean[0, :, :] - emis_JUN_mean[1, :, :]
bias_emis_JUL_mean = emis_JUL_mean[0, :, :] - emis_JUL_mean[1, :, :]



## carto ##
plt.ion()
plt.figure()
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., -30., 20))
m.drawmeridians(np.arange(-180., 180., 20))
xii,yii = m(*np.meshgrid(lon, lat))
clevs = (arange(-0.23, 0.18, 0.001))
cs = m.contourf(xii, yii, bias_emis_JUL_mean[:, :], clevs, cmap=cm.s3pcpn_l_r)
cbar = colorbar(cs)
cbar.set_label('ch1 - ch15')
plt.title('bias(emis) JUL - AMSUA')
xticks(arange(-180, 200, 20), rotation = 45)
yticks(arange(-90, -10, 20))


## selection des zones en fonction des valeurs de biais ##
bb0_JAN = nonzero((bias_emis_JAN_mean >= -0.01) & (bias_emis_JAN_mean <= 0.01))
bb1_JAN = nonzero((bias_emis_JAN_mean >= 0.10))
bb2_JAN = nonzero((bias_emis_JAN_mean <= -0.10))

plt.figure()
plt.plot(bias_emis_JAN_mean[ilat, :][bb0_JAN])

for ilon in range (0, len(lon)):
    ind = np.where((bias_emis_JAN_mean[:, ilon] >= -0.01) & (bias_emis_JAN_mean[:, ilon] <= 0.01))[0]      
    if len(ind) == 0:
        continue
    else: 
        a = ind
        print ilon, a

figure()