source: trunk/src/scripts_Laura/ARCTIC/Travail_CEN/map_monthly_mean_emis_parameters_whole_arctic.py @ 56

#!/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 arctic_map # function to regrid coast limits
import cartesian_grid_test # function to convert grid from polar to cartesian
import scipy.special
import ffgrid2
import map_ffgrid
from matplotlib import colors
from matplotlib.font_manager import FontProperties
import map_cartesian_grid


###############################
# time period characteristics #
###############################
MONTH = np.array(['01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12'])
month = np.array(['JANUARY', 'FEBRUARY', 'MARCH', 'APRIL', 'MAY', 'JUNE', 'JULY', 'AUGUST', 'SEPTEMBER', 'OCTOBER', 'NOVEMBER', 'DECEMBER'])
month_day = np.array([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31])
M = len(month)


########################
# grid characteristics #
########################
x0 = -3000. # min limit of grid
x1 = 2500. # max limit of grid
dx = 40.
xvec = np.arange(x0, x1+dx, dx)
nx = len(xvec) 
y0 = -3000. # min limit of grid
y1 = 3000. # max limit of grid
dy = 40.
yvec = np.arange(y0, y1+dy, dy)
ny = len(yvec)


grad_ratio = np.zeros([M, ny, nx], float)
std_grad_ratio = np.zeros([M, ny, nx], float)

spec_anom_23 = np.zeros([M, ny, nx], float)
std_spec_anom_23 = np.zeros([M, ny, nx], float)

spec_anom_89 = np.zeros([M, ny, nx], float)
std_spec_anom_89 = np.zeros([M, ny, nx], float)

ratio_anom_89 = np.zeros([M, ny, nx], float)
std_ratio_anom_89 = np.zeros([M, ny, nx], float)

std_spec_a = np.zeros([M, ny, nx], float)
std_spec_b = np.zeros([M, ny, nx], float)

for imo in range (0, M): 
    # AMSUA 89 std
    fichier_a = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_mean-std_data_lamb_spec_near_nadir_AMSUA89_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC')
    spec_a = fichier_a.variables['emis_spec_mean'][:]
    lamb_a = fichier_a.variables['emis_lamb_mean'][:]
    std_spec_a[imo, :, :] = fichier_a.variables['emis_spec_std'][:]
    #std_lamb_a[imo, :, :] = fichier_a.variables['emis_lamb_std'][:]
    fichier_a.close()
    # AMSUB 89 std
    fichier_b = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_mean-std_data_lamb_spec_near_nadir_AMSUB89_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC')
    spec_b = fichier_b.variables['emis_spec_mean'][:]
    lamb_b = fichier_b.variables['emis_lamb_mean'][:]
    std_spec_b[imo, :, :] = fichier_b.variables['emis_spec_std'][:]
    #std_lamb_b[imo, :, :] = fichier_b.variables['emis_lamb_std'][:]
    fichier_b.close()
    # read daily gradient ratio file
    print 'read daily gradient ratio for month ' + month[imo]
    fichier_grad_ratio = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_daily_grad_ratio_spec_23-89_near_nadir_AMSUA_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC')
    gr = fichier_grad_ratio.variables['grad_ratio'][:]
    fichier_grad_ratio.close()
    # read daily emis anomaly file for 23GHz
    print 'read daily emis anomaly 23GHz for month ' + month[imo]
    fichier_anom23 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_daily_data_lamb_spec_ratio_anomaly_near_nadir_AMSUA23_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC')
    sa_23 = fichier_anom23.variables['spec_anomaly'][:]
    fichier_anom23.close()
    # read daily emis anomaly file for 89GHz
    print 'read daily emis anomaly 89GHz for month ' + month[imo]
    fichier_anom89 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_daily_data_lamb_spec_ratio_anomaly_near_nadir_AMSUA89_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC')
    sa_89 = fichier_anom89.variables['spec_anomaly'][:]
    ra_89 = fichier_anom89.variables['ratio_anomaly'][:]
    fichier_anom89.close()
    print 'compute montly mean and std ' + month[imo]
    # calculate monthly means out of daily data
    for ilat in range (0, ny):
        for ilon in range (0, nx):
            # gradient ratio
            grad_ratio[imo, ilat, ilon] = mean(gr[0 : month_day[imo], ilat, ilon][nonzero(isnan(gr[0 : month_day[imo], ilat, ilon]) == False)])
            std_grad_ratio[imo, ilat, ilon] = sqrt((1. / (month_day[imo] - 1.)) * sum((gr[0 : month_day[imo], ilat, ilon][nonzero(isnan(gr[0 : month_day[imo], ilat, ilon]) == False)] - grad_ratio[imo, ilat, ilon])**2))
            # spec anomaly 23GHz
            spec_anom_23[imo, ilat, ilon] = mean(sa_23[0 : month_day[imo], ilat, ilon][nonzero(isnan(sa_23[0 : month_day[imo], ilat, ilon]) == False)])
            std_spec_anom_23[imo, ilat, ilon] = sqrt((1. / (month_day[imo] - 1.)) * sum((sa_23[0 : month_day[imo], ilat, ilon][nonzero(isnan(sa_23[0 : month_day[imo], ilat, ilon]) == False)] - spec_anom_23[imo, ilat, ilon])**2))
            # spec anomaly 89GHz
            spec_anom_89[imo, ilat, ilon] = mean(sa_89[0 : month_day[imo], ilat, ilon][nonzero(isnan(sa_89[0 : month_day[imo], ilat, ilon]) == False)])
            std_spec_anom_89[imo, ilat, ilon] = sqrt((1. / (month_day[imo] - 1.)) * sum((sa_89[0 : month_day[imo], ilat, ilon][nonzero(isnan(sa_89[0 : month_day[imo], ilat, ilon]) == False)] - spec_anom_89[imo, ilat, ilon])**2))
            # ratio anomaly 89GHz
            ratio_anom_89[imo, ilat, ilon] = mean(ra_89[0 : month_day[imo], ilat, ilon][nonzero(isnan(ra_89[0 : month_day[imo], ilat, ilon]) == False)])
            std_ratio_anom_89[imo, ilat, ilon] = sqrt((1. / (month_day[imo] - 1.)) * sum((ra_89[0 : month_day[imo], ilat, ilon][nonzero(isnan(ra_89[0 : month_day[imo], ilat, ilon]) == False)] - ratio_anom_89[imo, ilat, ilon])**2))
            # take out erroneous values of std (on land)
            if ((xvec[ilon] > -2400.) & (xvec[ilon] < -2160.) & (yvec[ilat] > 1240.) & (yvec[ilat] < 1560.)):
                grad_ratio[imo, ilat, ilon] = NaN
                spec_anom_23[imo, ilat, ilon] = NaN
                spec_anom_89[imo, ilat, ilon] = NaN
                ratio_anom_89[imo, ilat, ilon] = NaN
                std_spec_a[imo, ilat, ilon] = NaN
                std_spec_b[imo, ilat, ilon] = NaN


########################################################
# map correlation over the year between each parameter #
########################################################
corr_map1 = np.zeros([ny, nx], float)
corr_map2 = np.zeros([ny, nx], float)
corr_map3 = np.zeros([ny, nx], float)
corr_map4 = np.zeros([ny, nx], float)
corr_map5 = np.zeros([ny, nx], float)
corr_map6 = np.zeros([ny, nx], float)
corr_map7 = np.zeros([ny, nx], float)
corr_map8 = np.zeros([ny, nx], float)
corr_map9 = np.zeros([ny, nx], float)
corr_map10 = np.zeros([ny, nx], float)
for ilat in range (0, ny):
    for ilon in range (0, nx):
        corr_map1[ilat, ilon] = corrcoef(grad_ratio[:, ilat, ilon], spec_anom_23[:, ilat, ilon])[0][1]
        corr_map2[ilat, ilon] = corrcoef(grad_ratio[:, ilat, ilon], spec_anom_89[:, ilat, ilon])[0][1]
        corr_map3[ilat, ilon] = corrcoef(grad_ratio[:, ilat, ilon], ratio_anom_89[:, ilat, ilon])[0][1]
        corr_map4[ilat, ilon] = corrcoef(grad_ratio[:, ilat, ilon], std_spec_a[:, ilat, ilon])[0][1]
        corr_map5[ilat, ilon] = corrcoef(spec_anom_23[:, ilat, ilon], spec_anom_89[:, ilat, ilon])[0][1]
        corr_map6[ilat, ilon] = corrcoef(spec_anom_23[:, ilat, ilon], ratio_anom_89[:, ilat, ilon])[0][1]
        corr_map7[ilat, ilon] = corrcoef(spec_anom_23[:, ilat, ilon], std_spec_a[:, ilat, ilon])[0][1]
        corr_map8[ilat, ilon] = corrcoef(spec_anom_89[:, ilat, ilon], ratio_anom_89[:, ilat, ilon])[0][1]
        corr_map9[ilat, ilon] = corrcoef(spec_anom_89[:, ilat, ilon], std_spec_a[:, ilat, ilon])[0][1]
        corr_map10[ilat, ilon] = corrcoef(ratio_anom_89[:, ilat, ilon], std_spec_a[:, ilat, ilon])[0][1]


ion()
x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50()
x_coast = x_ind
y_coast = y_ind
z_coast = z_ind
map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, corr_map1[:, :], -1.1, 1.1, 0.1, cm.s3pcpn_l_r, 'correlation grad ratio - spec A23 anom')
plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA89_AMSUB89/maps/correlation_maps/correl_map_grad_ratio-param1.png')

map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, corr_map2[:, :], -1.1, 1.1, 0.1, cm.s3pcpn_l_r, 'correlation grad ratio - spec A89 anom')
plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA89_AMSUB89/maps/correlation_maps/correl_map_grad_ratio-param2.png')

map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, corr_map3[:, :], -1.1, 1.1, 0.1, cm.s3pcpn_l_r, 'correlation grad ratio - ratio A89 anom')
plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA89_AMSUB89/maps/correlation_maps/correl_map_grad_ratio-param3.png')

map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, corr_map4[:, :], -1.1, 1.1, 0.1, cm.s3pcpn_l_r, 'correlation grad ratio - spec std A89')
plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA89_AMSUB89/maps/correlation_maps/correl_map_grad_ratio-param4.png')

map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, corr_map5[:, :], -1.1, 1.1, 0.1, cm.s3pcpn_l_r, 'correlation spec A23 anom - spec A89 anom')
plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA89_AMSUB89/maps/correlation_maps/correl_map_grad_ratio-param5.png')

map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, corr_map6[:, :], -1.1, 1.1, 0.1, cm.s3pcpn_l_r, 'correlation spec A23 anom - ratio A89 anom')
plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA89_AMSUB89/maps/correlation_maps/correl_map_grad_ratio-param6.png')

map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, corr_map7[:, :], -1.1, 1.1, 0.1, cm.s3pcpn_l_r, 'correlation spec A23 anom - spec std A89')
plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA89_AMSUB89/maps/correlation_maps/correl_map_grad_ratio-param7.png')

map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, corr_map8[:, :], -1.1, 1.1, 0.1, cm.s3pcpn_l_r, 'correlation spec A89 anom - ratio A89 anom')
plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA89_AMSUB89/maps/correlation_maps/correl_map_grad_ratio-param8.png')

map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, corr_map9[:, :], -1.1, 1.1, 0.1, cm.s3pcpn_l_r, 'correlation spec A89 anom - std spec A89')
plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA89_AMSUB89/maps/correlation_maps/correl_map_grad_ratio-param9.png')

map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, corr_map10[:, :], -1.1, 1.1, 0.1, cm.s3pcpn_l_r, 'correlation ratio A89 anom - std spec A89')
plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA89_AMSUB89/maps/correlation_maps/correl_map_grad_ratio-param10.png')


'''
# test
map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, std_spec_a[imo, :, :], 0., 0.12, 0.001, cm.s3pcpn_l_r, 'test')
'''

############################################################
# correlation matrix between each parameter for each month #
############################################################
# reshape matrix into vector fo calculation or correlation per month between each parameter
corr_mat = np.zeros([M, 5, 5], float)
for imo in range (0, M):
    print 'month ' + month[imo]
    a = reshape(std_spec_a[imo, :, :], size(std_spec_a[imo, :, :]))
    b = reshape(std_spec_b[imo, :, :], size(std_spec_b[imo, :, :]))
    c = reshape(grad_ratio[imo, :, :], size(grad_ratio[imo, :, :]))
    d = reshape(spec_anom_23[imo, :, :], size(spec_anom_23[imo, :, :]))
    e = reshape(spec_anom_89[imo, :, :], size(spec_anom_89[imo, :, :]))
    f = reshape(ratio_anom_89[imo, :, :], size(ratio_anom_89[imo, :, :]))
    aa = a[nonzero(isnan(a) == False)]
    bb = b[nonzero(isnan(b) == False)]
    cc = c[nonzero(isnan(c) == False)]
    dd = d[nonzero(isnan(d) == False)]
    ee = e[nonzero(isnan(e) == False)]
    ff = f[nonzero(isnan(f) == False)]
    print len(aa), len(bb), len(cc), len(dd), len(ee), len(ff)
    params = np.array([aa, cc, dd, ee, ff])
    for ii in range (0, 5):
        for jj in range (0, 5):
            corr_mat[imo, ii, jj] = corrcoef(params[ii, :], params[jj, :])[0][1]


ion()
for imo in range (0, M):
    figure()
    pc = pcolor(corr_mat[imo, :, :], vmin = -1., vmax = 1., cmap = 'RdBu')
    cbar = colorbar(pc)
    xticks(np.arange(0.5, 5.5, 1.), np.array(['std spec A89', 'grad ratio A', 'spec A23 anom', 'spec A89 anom', 'ratio A89 anom']), rotation = 20)
    yticks(np.arange(0.5, 5.5, 1.), np.array(['std spec A89', 'grad ratio A', 'spec A23 anom', 'spec A89 anom', 'ratio A89 anom']), rotation = 70)
    cbar.set_label('correlation')
    title(month[imo] + ' 2009')
    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA89_AMSUB89/correlation_matrix/corr_matrix_emis_params_temporal_std89_' + MONTH[imo] + month[imo] + '2009.png')


'''
############################
# map monthly mean and std #
############################
#ion()
x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50()
x_coast = x_ind
y_coast = y_ind
z_coast = z_ind
for imo in range (0, M):
    print 'map month ' + month[imo]
    # emis spec std AMSUA89
    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, std_spec_a[imo, :, :], 0., 0.12, 0.001, cm.s3pcpn_l_r, 'emis spec monthly std')
    title('AMSUA 89GHz - ' + month[imo] + ' 2009')
    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/comparison_lamb_spec/space_evolution/EMIS/cartesian_grid/monthly_std/std_emis_spec_AMSUA89_' + MONTH[imo] + month[imo] + '2009.png')
    # emis spec std AMSUB89
    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, std_spec_b[imo, :, :], 0., 0.12, 0.001, cm.s3pcpn_l_r, 'emis spec monthly std')
    title('AMSUB 89GHz - ' + month[imo] + ' 2009')
    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/comparison_lamb_spec/space_evolution/EMIS/cartesian_grid/monthly_std/std_emis_spec_AMSUB89_' + MONTH[imo] + month[imo] + '2009.png')
    # gradient ratio
    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, grad_ratio[imo, :, :], grad_ratio[nonzero(isnan(grad_ratio) == False)].min(), grad_ratio[nonzero(isnan(grad_ratio) == False)].max(), 0.001, cm.s3pcpn_l_r, 'Gradient ratio monthly mean')
    title('AMSUA - ' + month[imo] + ' 2009')
    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/maps/grad_ratio/map_monthly_mean_grad_ratio_AMSUA23-89_' + month[imo] + '2009.png')
    # spec anomaly 23GHz
    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, spec_anom_23[imo, :, :], spec_anom_23[nonzero(isnan(spec_anom_23) == False)].min(), spec_anom_23[nonzero(isnan(spec_anom_23) == False)].max(), 0.001, cm.s3pcpn_l_r, 'Emis spec anomaly 23GHz monthly mean')
    title('AMSUA - ' + month[imo] + ' 2009')
    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/maps/spec_anom_23/map_monthly_mean_spec_anomaly_AMSUA23_' + month[imo] + '2009.png')
    # spec anomaly 89GHz
    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, spec_anom_89[imo, :, :], spec_anom_89[nonzero(isnan(spec_anom_89) == False)].min(), spec_anom_89[nonzero(isnan(spec_anom_89) == False)].max(), 0.001, cm.s3pcpn_l_r, 'Emis spec anomaly 89GHz monthly mean')
    title('AMSUA - ' + month[imo] + ' 2009')
    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/maps/spec_anom_89/map_monthly_mean_spec_anomaly_AMSUA89_' + month[imo] + '2009.png')
    # ratio anomaly 89GHz
    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, ratio_anom_89[imo, :, :], ratio_anom_89[nonzero(isnan(ratio_anom_89) == False)].min(), ratio_anom_89[nonzero(isnan(ratio_anom_89) == False)].max(), 0.001, cm.s3pcpn_l_r, 'Emis ratio anomaly 89GHz monthly mean')
    title('AMSUA - ' + month[imo] + ' 2009')
    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/maps/ratio_anom_89/map_monthly_mean_ratio_anomaly_AMSUA89_' + month[imo] + '2009.png')
'''


'''
map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, std_spec_a[imo, :, :], std_spec_a[imo, :, :][nonzero(std_spec_a[imo, :, :] != 0.)].min(), std_spec_a[imo, :, :][nonzero(std_spec_a[imo, :, :] != 0.)].max(), 0.001, cm.s3pcpn_l_r, 'Gradient ratio monthly std')

map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, gr[15, :, :], gr[15, :, :][nonzero(isnan(gr[15, :, :]) == False)].min(), gr[15, :, :][nonzero(isnan(gr[15, :, :]) == False)].max(), 0.001, cm.s3pcpn_l_r, 'daily gradient ratio (01-01-2009)')
contour(xvec[xxi:xxf+1], yvec[yyi:yyf+1], )
'''