#!/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)




##################################################################
# compute monthly means of emissivity parameters on sea ice area #
##################################################################
spec_emis89 = np.zeros([M, ny, nx], float)
spec_emis23 = np.zeros([M, ny, nx], float)

#std_spec_emis89 = np.zeros([M, ny, nx], float)
#std_spec_emis23 = np.zeros([M, ny, nx], float)

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

spec_anom89 = np.zeros([M, ny, nx], float)
ratio_anom89 = np.zeros([M, ny, nx], float)

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

for imo in range (0, M): 
    print 'compute emissivity parameters for month ' + month[imo]
    # open daily spec emissivity at 89GHz
    fichier_daily89 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/sub_classification/cartesian_grid_map_sea_ice_extent_with-AMSUA23-and-89_' + month[imo] + '2009_AMSUA89_spec_thresholds.nc', 'r', format='NETCDF3_CLASSIC')
    a = fichier_daily89.variables['spec_ice_area'][:]
    fichier_daily89.close()
    # open daily spec emissivity at 23GHz
    fichier_daily23 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/sub_classification/cartesian_grid_map_sea_ice_extent_with-AMSUA23-and-89_' + month[imo] + '2009_AMSUA23_spec_thresholds.nc', 'r', format='NETCDF3_CLASSIC')
    b = fichier_daily23.variables['spec_ice_area'][:]
    fichier_daily23.close()
    ## compute monthly mean of spec emissivity at 89GHz
    print 'monthly spec emis'
    for ilon in range (0, nx):
        for ilat in range (0, ny):
            spec_emis89[imo, ilat, ilon] = mean(a[ilat, ilon, 0 : month_day[imo]][nonzero(isnan(a[ilat, ilon, 0 : month_day[imo]]) == False)])
            #std_spec_emis89[imo, ilat, ilon] = sqrt((1. / (month_day[imo] - 1.)) * sum((a[ilat, ilon, 0 : month_day[imo]][nonzero(isnan(a[ilat, ilon, 0 : month_day[imo]]) == False)] - spec_emis89[imo, ilat, ilon])**2.))
            spec_emis23[imo, ilat, ilon] = mean(b[ilat, ilon, 0 : month_day[imo]][nonzero(isnan(b[ilat, ilon, 0 : month_day[imo]]) == False)])
            #std_spec_emis23[imo, ilat, ilon] = sqrt((1. / (month_day[imo] - 1.)) * sum((b[ilat, ilon, 0 : month_day[imo]][nonzero(isnan(b[ilat, ilon, 0 : month_day[imo]]) == False)] - spec_emis23[imo, ilat, ilon])**2.))
    # compute monthly gradient ratio 23-89GHz
    print 'monthly gradient ratio 23-89GHz'
    grad_ratio[imo, :, :] = spec_emis23[imo, :, :] - spec_emis89[imo, :, :]
    # open anomaly of spec emissivity and emissivity ratio anomaly at 89GHz
    fichier_anom89 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/sub_classification/cartesian_grid_map_sea_ice_extent_with-AMSUA23-and-89_param_anomaly_' + month[imo] + '2009_AMSUA89.nc', 'r', format = 'NETCDF3_CLASSIC')
    spec_anom89[imo, :, :] = fichier_anom89.variables['spec_anomaly'][:]
    ratio_anom89[imo, :, :] = fichier_anom89.variables['ratio_anomaly'][:]
    fichier_anom89.close()
    # open anomaly of spec emissivity anomaly at 23GHz
    fichier_anom23 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/sub_classification/cartesian_grid_map_sea_ice_extent_with-AMSUA23-and-89_param_anomaly_' + month[imo] + '2009_AMSUA23.nc', 'r', format = 'NETCDF3_CLASSIC')
    spec_anom23[imo, :, :] = fichier_anom23.variables['spec_anomaly'][:]
    fichier_anom23.close()
    


cumul = abs(spec_anom89) + abs(spec_anom23) + abs(ratio_anom89) + abs(grad_ratio)
cumul_index = np.zeros([M, ny, nx], float)
for imo in range (0, M):
    print 'month ' + month[imo]
    max_cumul = cumul[imo, :, :][nonzero(isnan(cumul[imo, :, :]) == False)].max()
    for ilon in range (0, nx):
        for ilat in range (0, ny):
            cumul_index[imo, ilat, ilon] = cumul[imo, ilat, ilon] / max_cumul


    


##############################################
# maps monthly mean of emissivity parameters #
##############################################
#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]
    # cumulation index
    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, cumul_index[imo, :, :], 0., 1., 0.01, cm.s3pcpn_l_r, 'monthly cumulation index')
    title('AMSUA - ' + month[imo] + ' 2009')
    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/maps/cumul_params/cumul_all_parameters/map_cumulation_index_' + MONTH[imo] + month[imo] + '2009_sea_ice_extent.png')
    '''
    # emis spec std 23GHz
    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, std_spec_emis23[imo, :, :], 0., 0.12, 0.001, cm.s3pcpn_l_r, 'monthly emis spec std')
    title('AMSUA 23GHz - ' + month[imo] + ' 2009')
    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/maps/emiss_spec/std/std_emis_spec_map_AMSUA23_' + MONTH[imo] + month[imo] + '_2009_arctic_sea_ice.png')
    # emis spec std 89GHz
    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, std_spec_emis89[imo, :, :], 0., 0.2, 0.001, cm.s3pcpn_l_r, 'monthly emis spec std')
    title('AMSUA 89GHz - ' + month[imo] + ' 2009')
    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/maps/emiss_spec/std/std_emis_spec_map_AMSUA89_' + MONTH[imo] + month[imo] + '_2009_arctic_sea_ice.png')
    # emis spec anomaly AMSUA23
    print 'map emis spec anomaly AMSUA23'
    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, spec_anom23[imo, :, :], -0.1, 0.2, 0.01, cm.s3pcpn_l_r, 'monthly emis spec anomaly')
    title('AMSUA 23GHz - ' + month[imo] + ' 2009')
    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/maps/emiss_anomaly/spec/anomaly_spec_map_AMSUA23_' + MONTH[imo] + month[imo] + '_2009_arctic_sea_ice.png')
    # emis spec anomaly AMSUA89
    print 'map emis spec anomaly AMSUA89'
    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, spec_anom89[imo, :, :], -0.2, 0.2, 0.01, cm.s3pcpn_l_r, 'monthly emis spec anomaly')
    title('AMSUA 89GHz - ' + month[imo] + ' 2009')
    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/maps/emiss_anomaly/spec/anomaly_spec_map_AMSUA89_' + MONTH[imo] + month[imo] + '_2009_arctic_sea_ice.png')
    # emis gradient ratio
    print 'map gradient ratio'
    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, grad_ratio[imo, :, :], -0.2, 0.3, 0.01, cm.s3pcpn_l_r, 'monthly gradient ratio')
    title('AMSUA [23-89]GHz - ' + month[imo] + ' 2009')
    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/maps/gradient_ratio_23-89/gradient_ratio_map_AMSUA23-89_spec_' + MONTH[imo] + month[imo] + '_2009_arctic_sea_ice.png')
    # emis ratio anomaly AMSUA89
    print 'map emis ratio anomaly AMSUA89'
    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, ratio_anom89[imo, :, :], -2, 2, 0.1, cm.s3pcpn_l_r, 'monthly emis ratio anomaly')
    title('AMSUA 89GHz - ' + month[imo] + ' 2009')
    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/maps/emiss_anomaly/ratio/anomaly_ratio_map_AMSUA89_' + MONTH[imo] + month[imo] + '_2009_arctic_sea_ice.png')
    '''