source: trunk/src/scripts_Laura/ARCTIC/Travail_CEN/daily_param_emissivity_all_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)




###################################

frequ = np.array(['23', '89'])




emis_spec_anom = np.zeros([len(frequ), M, 31, ny, nx], float)
emis_lamb_anom = np.zeros([len(frequ), M, 31, ny, nx], float)
emis_diff_anom = np.zeros([len(frequ), M, 31, ny, nx], float)
emis_ratio_anom = np.zeros([len(frequ), M, 31, ny, nx], float)
for ifr in range (1, len(frequ)):
    print 'frequency ' + frequ[ifr]
    for imo in range (0, M):
        print 'month ' + month[imo]
        print 'read daily emis spec lamb and diff'
        fichier_glob = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA' + str(frequ[ifr]) + '_' + month[imo] + '2009.nc', 'r', format = 'NETCDF3_CLASSIC')
        emis_spec = fichier_glob.variables['e_spec'][:]
        emis_lamb = fichier_glob.variables['e_lamb'][:]
        emis_diff = fichier_glob.variables['e_spec_lamb'][:]
        fichier_glob.close()
        print 'calculate daily emis_ratio'
        emis_ratio = np.zeros([ny, nx, month_day[imo]], float)
        for ijr in range (0, month_day[imo]):
            for ilat in range (0, ny):
                for ilon in range (0, nx):
                    emis_ratio[ilat, ilon, ijr] = ((emis_lamb[ilat, ilon, ijr] - emis_spec[ilat, ilon, ijr]) / emis_spec[ilat, ilon, ijr]) * 100.
        print 'calculate daily anomaly'
        a = np.zeros([month_day[imo], ny * nx], float)
        b = np.zeros([month_day[imo], ny * nx], float)
        c = np.zeros([month_day[imo], ny * nx], float)
        d = np.zeros([month_day[imo], ny * nx], float)
        for ijr in range (0, month_day[imo]):
            print 'date ' + str(ijr + 1) + ' ' + month[imo] + ' 2009'
            a[ijr, :] = reshape(emis_spec[:, :, ijr], size(emis_spec[:, :, ijr]))
            b[ijr, :] = reshape(emis_lamb[:, :, ijr], size(emis_lamb[:, :, ijr]))
            c[ijr, :] = reshape(emis_diff[:, :, ijr], size(emis_diff[:, :, ijr]))
            d[ijr, :] = reshape(emis_ratio[:, :, ijr], size(emis_ratio[:, :, ijr]))
            for ilat in range (0, ny):
                for ilon in range (0, nx):
                    emis_spec_anom[ifr, imo, ijr, ilat, ilon] = emis_spec[ilat, ilon, ijr] - mean(a[ijr, :][nonzero(isnan(a[ijr, :]) == False)])
                    emis_lamb_anom[ifr, imo, ijr, ilat, ilon] = emis_lamb[ilat, ilon, ijr] - mean(b[ijr, :][nonzero(isnan(b[ijr, :]) == False)])
                    emis_diff_anom[ifr, imo, ijr, ilat, ilon] = emis_diff[ilat, ilon, ijr] - mean(c[ijr, :][nonzero(isnan(c[ijr, :]) == False)])
                    emis_ratio_anom[ifr, imo, ijr, ilat, ilon] = emis_ratio[ilat, ilon, ijr] - mean(d[ijr, :][nonzero(isnan(d[ijr, :]) == False)])
        ########################
        # stack in netcdf file #
        ######################## 
        print 'stack data in file - month ' + str(month[imo])
        rootgrp = 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_AMSUA' + str(frequ[ifr]) + '_' + month[imo] + '2009.nc', 'w', format='NETCDF3_CLASSIC')
        rootgrp.createDimension('longitude', nx)
        rootgrp.createDimension('latitude', ny)
        rootgrp.createDimension('days', month_day[imo])
        nc_lon = rootgrp.createVariable('longitude', 'f', ('longitude',))
        nc_lat = rootgrp.createVariable('latitude', 'f', ('latitude',))
        nc_days = rootgrp.createVariable('days', 'f', ('days',))
        nc_anom_spec = rootgrp.createVariable('spec_anomaly', 'f', ('days', 'latitude', 'longitude'))
        nc_anom_lamb = rootgrp.createVariable('lamb_anomay', 'f', ('days', 'latitude', 'longitude'))
        nc_anom_diff = rootgrp.createVariable('diff_anomaly', 'f', ('days', 'latitude', 'longitude'))
        nc_anom_ratio = rootgrp.createVariable('ratio_anomaly', 'f', ('days', 'latitude', 'longitude'))
        nc_lon[:] = xvec
        nc_lat[:] = yvec
        nc_days[:] = np.arange(0, month_day[imo], 1)
        nc_anom_spec[:] = emis_spec_anom[ifr, imo, 0 : month_day[imo], :, :]
        nc_anom_lamb[:] = emis_lamb_anom[ifr, imo,0 : month_day[imo], :, :]
        nc_anom_diff[:] = emis_diff_anom[ifr, imo,0 : month_day[imo], :, :]
        nc_anom_ratio[:] = emis_ratio_anom[ifr, imo,0 : month_day[imo], :, :]
        rootgrp.close()