source: trunk/src/scripts_Laura/ARCTIC/Travail_CEN/daily_evolution_zones_osisaf_emis_spec_lamb.py @ 55

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



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 parameters #
###################
dx=1.
dy=0.5
x0, x1 = -180., 180.
y0, y1 = 30., 90.
xvec = np.arange(x0, x1 + dx, dx)
yvec = np.arange(y0, y1 + dy, dy)
nx = len(xvec)
ny = len(yvec)



########################################
# different ice types from OSISAF data #
########################################
ICE = np.array(['Open_Sea', 'Seasonal_Ice', 'Permanent_Ice', 'Ambiguous_Ice_Type'])



############################
# definition of study zone #
############################
## latitude ##
zlat = np.where((yvec >= 80.) & (yvec <= 83.))[0]
zlat_vec = yvec[zlat]
nlat = len(zlat)
i0 = zlat[0]
i1 = zlat[nlat - 1] + 1
## longitude ##
zlon = np.where((xvec <= -100.) & (xvec >= -130.))[0]
zlon_vec = xvec[zlon]
nlon = len(zlon)
j0 = zlon[0]
j1 = zlon[nlon - 1] + 1


##########################################
# calculation of mean SPEC and LAMB emis #
##########################################
type_map = np.zeros([M, ny, nx], float)
mean_spec_m = np.zeros([M, 31], float)
ect_spec_m = np.zeros([M, 31], float)
mean_lamb_m = np.zeros([M, 31], float)
ect_lamb_m = np.zeros([M, 31], float)
pred_oit_m = np.zeros([M, 31], float)
ii = np.zeros([M], float)
for imo in range (0, M):
    print 'month' + month[imo]
    ############################################
    # Read daily gridded OSISAF and AMSUB data #
    ############################################
    fichier = Dataset('/mma/hermozol/Documents/Data/daily_OSISAF_SPEC_LAMB/daily_OSISAF_emis_spec_lamb_AMSUB_ffgrid_' + str(month[imo]) + '_2009.nc', 'r', format='NETCDF3_CLASSIC')
    lon = fichier.variables['longitude'][:]
    lat = fichier.variables['latitude'][:]
    day = fichier.variables['days'][:]
    osi_type = fichier.variables['osi_ice_type'][:]
    type_map[imo, :, :] = osi_type[0, :, :]
    amsub_spec = fichier.variables['spec_emis'][:]
    amsub_lamb = fichier.variables['lamb_emis'][:]
    fichier.close()
    mean_spec = np.zeros([month_day[imo]], float)
    ect_spec = np.zeros([month_day[imo]], float)
    mean_lamb = np.zeros([month_day[imo]], float)
    ect_lamb = np.zeros([month_day[imo]], float)
    pred_oit = np.zeros([month_day[imo]], float)
    for ijr in range (0, month_day[imo]):
        s = reshape(amsub_spec[ijr, i0 : i1, j0 : j1], size(amsub_spec[ijr, i0 : i1, j0 : j1]))
        l = reshape(amsub_lamb[ijr, i0 : i1, j0 : j1], size(amsub_lamb[ijr, i0 : i1, j0 : j1]))
        oit = reshape(osi_type[ijr, i0 : i1, j0 : j1], size(osi_type[ijr, i0 : i1, j0 : j1]))
        ## take out emissivities == 0.
        for ks in range (0, len(s)):
            if (s[ks] == 0.):
                s[ks] = NaN
        ## take out emissivities == 0.
        for kl in range (0, len(l)):
            if (l[kl] == 0.):
                l[kl] = NaN
        mean_spec[ijr] = mean(s[nonzero(isnan(s) == False)])
        ect_spec[ijr] = (1. / (nlon * nlat - 1)) * sqrt(sum((s[nonzero(isnan(s) == False)] - mean_spec[ijr]) ** 2))
        if (ect_spec[ijr] == 0.):
            ect_spec[ijr] = NaN
        mean_lamb[ijr] = mean(l[nonzero(isnan(l) == False)])
        ect_lamb[ijr] = (1. / (nlon * nlat - 1)) * sqrt(sum((l[nonzero(isnan(l) == False)] - mean_lamb[ijr]) ** 2))
        if (ect_lamb[ijr] == 0.):
            ect_lamb[ijr] = NaN         
        print 'day' + str(ijr + 1)
        pred_oit[ijr] = round(mean(oit[nonzero(isnan(oit) == False)])) - 1
    mean_spec_m[imo, 0 : month_day[imo]] = mean_spec[:]
    ect_spec_m[imo, 0 : month_day[imo]] = ect_spec[:]
    mean_lamb_m[imo, 0 : month_day[imo]] = mean_lamb[:]
    ect_lamb_m[imo, 0 : month_day[imo]] = ect_lamb[:]
    ii[imo] = int(pred_oit[0])


## map of study zone ##
for imo in range (0, M):
    map_ffgrid.draw_map_l (xvec, yvec, 0, 5, 1, type_map[8, :, :], colors.ListedColormap(['1.', '0.25', '0.50', '0.75']), 'OSISAF Ice Type')
    title(month[imo] + '2009')






SPEC = np.zeros([M], float)
LAMB = np.zeros([M], float)
ECT_SPEC = np.zeros([M], float)
ECT_LAMB = np.zeros([M], float)
for imo in range (0, M):
    #map_ffgrid.draw_map_l (xvec, yvec, 0, 5, 1, type_map[1, :, :], colors.ListedColormap(['1.', '0.25', '0.50', '0.75']), 'OSISAF Ice Type')
    SPEC[imo] = mean(mean_spec_m[imo, :][nonzero(isnan(mean_spec_m[imo, :]) == False)])
    ECT_SPEC[imo] = mean(ect_spec_m[imo, :][nonzero(isnan(ect_spec_m[imo, :]) == False)])
    LAMB[imo] = mean(mean_lamb_m[imo, :][nonzero(isnan(mean_lamb_m[imo, :]) == False)])
    ECT_LAMB[imo] = mean(ect_lamb_m[imo, :][nonzero(isnan(ect_lamb_m[imo, :]) == False)])



plt.ion()
figure()
errorbar(np.arange(0, M, 1), SPEC, label = 'SPEC', yerr = ECT_SPEC * 5)
errorbar(np.arange(0, M, 1), LAMB, label = 'LAMB', yerr = ECT_LAMB * 5)
legend(loc = 3)
ylabel('emissivity and std (x5)')
xlim(-1, M)
ylim(0.5, 1.)
xticks(np.arange(0, M, 1), month, rotation = 25)
grid(True)
twinx()
plot(ii, '--k', label = 'OSISAF ice type')
legend(loc = 1)
yticks(np.arange(0, 4, 1), np.array(['ice free', 'seasonal', 'permanent', 'amiguous']), rotation = 60)
ylim(0 , 4)
title('zone : OSISAF 2009 - AMSUB')
plt.savefig('/mma/hermozol/Documents/figure_output/comp_lamb_spec/time_evolution/daily_evo_emis_spec_lamb_' + ICE[ii[0]] + '_' + month[imo] + '2009_AMSUB.png')

















#########
# tests #
#########
z0 = spec[0, i0 : i1, j0 : j1].min()
z1 = spec[0, i0 : i1, j0 : j1].max()
map_ffgrid.draw_map_l (zlon_vec, zlat_vec, 0.4, 1.02, 0.01, spec[0, i0 : i1, j0 : j1] , cm.s3pcpn_l_r, 'emis SPEC')
cmap2 = colors.ListedColormap(['1.', '0.25', '0.50', '0.75'])
map_ffgrid.draw_map_l (zlon_vec, zlat_vec, 0, 5, 1, osi_type[6, 0, i0 : i1, j0 : j1] , cmap2, 'OSISAF ice type')


map_ffgrid.draw_map_l (xvec, yvec, 0.5, 1.02, 0.01, spec[6, 0, :, :], cm.s3pcpn_l_r, 'emissivity SPEC')
map_ffgrid.draw_map_l (xvec, yvec, 0.5, 1.02, 0.01, lamb[6, 0, :, :], cm.s3pcpn_l_r, 'emissivity LAMB')
map_ffgrid.draw_map_l (xvec, yvec, 0, 5, 1, type_map[1, :, :], colors.ListedColormap(['1.', '0.25', '0.50', '0.75']), 'OSISAF Ice Type')
cbar.set_ticklabels(['no ice','seasonal ice','permanent ice','ambiguous ice type'])
#xii, yii = m(*np.meshgrid(xvec, yvec))
#cmap2 = colors.ListedColormap(['1.', '0.25', '0.50', '0.75'])
#m.contour(xii, yii, osi_ice_type[0, :, :], np.arange(0, 5, 1))
#CS = plt.contour(xii, yii, osi_ice_type[0, :, :], 3,
                 colors='k', # negative contours will be dashed by default
                 )