source: trunk/src/scripts_Laura/ARCTIC/Travail_CEN/discrim_rate_espec-elamb_espec.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
import map_cartesian_grid

    

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)



############################################
# Read daily gridded OSISAF and AMSUB data #
############################################
sl_rate = np.zeros([M, 31, ny, nx], float)
for imo in range (0, M):
    print 'month: ' + month[imo]
    ### AMSUA23 ###
    print 'open data file'
    fichier = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA30_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC')
    xdist = fichier.variables['longitude'][:]
    ydist = fichier.variables['latitude'][:]
    day = fichier.variables['days'][:]
    emis_spec = fichier.variables['e_spec'][:]
    emis_lamb = fichier.variables['e_lamb'][:]
    fichier.close()
    print 'calculate rate'
    for ijr in range (0, month_day[imo]):
        for ilat in range (0, len(ydist)):
            for ilon in range (0, len(xdist)):
                sl_rate[imo, ijr, ilat, ilon] = ((emis_lamb[ilat, ilon, ijr] - emis_spec[ilat, ilon, ijr]) / emis_spec[ilat, ilon, ijr]) * 100.
    



#############################################
# compute monthly means of emissivity ratio #
#############################################
print 'compute monthly mean of emissivity ratio'
monthly_ratio = np.zeros([M, ny, nx], float)
for imo in range (0, M):
    print 'month ' + str(month[imo])
    for ilon in range (0, len(xdist)):
        for ilat in range (0, len(ydist)):
            monthly_ratio[imo, ilat, ilon] = mean(sl_rate[imo, 0 : month_day[imo], ilat, ilon][nonzero(isnan(sl_rate[imo, 0 : month_day[imo], ilat, ilon]) == False)])
    

    

    

#########################################
# map monthly means of emissivity ratio #
#########################################
print 'start mapping'
## parameters of coast coordinates for mapping ##
plt.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
colormap = cm.s3pcpn_l_r
## start mapping ratio ##
for imo in range (0, M):
    print 'map month ' +str(month[imo])
    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xdist, ydist, monthly_ratio[imo, :, :], 0., 6., 0.1, colormap, 'emissivity ratio [(LAMB - SPEC)/SPEC]*100')
    title(month[imo] + ' 2009 - AMSUA 30GHz')
    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/comparison_lamb_spec/space_evolution/EMIS/cartesian_grid/emis_ratio/emis_ratio_spec-lamb_AMSUA30_' + month[imo] + '2009.png')






###########################################################
# stack monthly means of emissivity ratio in netcdf files #
###########################################################
print 'start stacking'
## AMSUA ratio ##
for imo in range (0, M):
    print 'stack in file month ' + str(month[imo])
    rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_lamb-spec_ratio_near_nadir_AMSUA30_' + month[imo] + '2009.nc', 'w', format='NETCDF3_CLASSIC')
    rootgrp.createDimension('longitude', len(xdist))
    rootgrp.createDimension('latitude', len(ydist))
    nc_lon = rootgrp.createVariable('longitude', 'f', ('longitude',))
    nc_lat = rootgrp.createVariable('latitude', 'f', ('latitude',))
    nc_ratio = rootgrp.createVariable('emis_ratio', 'f', ('latitude', 'longitude'))
    nc_lon[:] = xdist
    nc_lat[:] = ydist
    nc_ratio[:] = monthly_ratio[imo, :, :]
    rootgrp.close()