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




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)

########################### for AMSUA (4 DIFFERENT FREQUENCIES) #########################

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

########################################
# reads NETCFD files of sea ice extent #
########################################
ice_spec = np.zeros([4, M, ny, nx], float)
ice_lamb = np.zeros([4, M, ny, nx], float)

vec_read_month = np.arange(0, 600, 50)

spec_lim = np.zeros([4, M], float)
lamb_lim = np.zeros([4, M], float)

for ifr in range (0, 4):
    print 'frequency ' + frequ[ifr]
    fichierA_dat = open('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/AMSUA' + frequ[ifr] + '_data_classification_parameters_ice_no-ice_2009.dat','r')
    numlines = 0
    for line in fichierA_dat: numlines += 1
    fichierA_dat.close()
    fichierA_dat = open('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/AMSUA' + frequ[ifr] + '_data_classification_parameters_ice_no-ice_2009.dat','r')
    nbtotal=numlines-1
    iligne=0
    mo = np.zeros([nbtotal],object) # month
    spec_emis_thresh = np.zeros([nbtotal],float) # emissivity threshold with spec emis
    lamb_emis_thresh = np.zeros([nbtotal],float) # emissivity threshold with spec emis
    while (iligne < nbtotal) :
         line = fichierA_dat.readline()
         liste = line.split()
         mo[iligne] = str(liste[0])
         spec_emis_thresh[iligne] = float(liste[7])
         lamb_emis_thresh[iligne] = float(liste[8])
         iligne=iligne+1
    fichierA_dat.close()
    spec_lim[ifr, :] = spec_emis_thresh[vec_read_month] # read spec emis threshold for each month and each frequence
    lamb_lim[ifr, :] = lamb_emis_thresh[vec_read_month] # read lamb emis threshold for each month and each frequence
    for imo in range (0, M):
        print 'month ' + month[imo]
        # AMSUA NETCDF file
        fichierA_nc = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/cartesian_grid_map_ice_no-ice_' + month[imo] + '2009_AMSUA' + frequ[ifr] + '_spec_lamb_thresholds.nc', 'r', format='NETCDF3_CLASSIC')
        ice_spec[ifr, imo, :, :] = fichierA_nc.variables['spec_ice_area'][:]
        ice_lamb[ifr, imo, :, :] = fichierA_nc.variables['lamb_ice_area'][:]
        fichierA_nc.close()
        



##############################################################
# map sea ice extent with lamb and spec emissivity thresholds #
##############################################################    
print 'start mapping'
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 ifr in range (0, 4):
    print 'map for frequency ' + frequ[ifr] + 'GHZ'
    for imo in range (0, M):
        # spec threshold
        print 'map with spec threshold month ' + month[imo]
        map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, new_ice_spec[ifr, imo, :, :], 0, 2, 1, colors.ListedColormap(['blue']), 'Sea ice extent with emis spec > ' + str(spec_lim[ifr, imo])[0:5])
        title('AMSUA ' + frequ[ifr] + 'GHz - ' + month[imo] + ' 2009')
        plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/maps/sea_ice_extent/corrected_extent/AMSUA' + frequ[ifr] + '_ice_emis_spec_thresh_' + MONTH[imo] + month[imo] + '2009.png')
        # lamb threshold
        print 'map with lamb threshold month ' + month[imo]
        map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, ice_lamb[ifr, imo, :, :], -1, 2, 1, colors.ListedColormap(['0.5', 'blue']), 'Sea ice extent with emis lamb > ' + str(lamb_lim[ifr, imo])[0:5])
        title('AMSUA ' + frequ[ifr] + 'GHz - ' + month[imo] + ' 2009')
        plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/maps/sea_ice_extent/AMSUA' + frequ[ifr] + '_ice_emis_lamb_thresh_' + MONTH[imo] + month[imo] + '2009.png')
        




########################### for AMSUB (1 FREQUENCY 89GHz) #########################


########################################
# reads NETCFD files of sea ice extent #
########################################
ice_spec = np.zeros([M, ny, nx], float)
ice_lamb = np.zeros([M, ny, nx], float)

vec_read_month = np.arange(0, 600, 50)

spec_lim = np.zeros([M], float)
lamb_lim = np.zeros([M], float)

fichierB_dat = open('/net/argos/data/parvati/lahlod/ARCTIC/AMSUB_ice_class/AMSUB89_data_classification_parameters_ice_no-ice_2009.dat','r')
numlines = 0
for line in fichierB_dat: numlines += 1

fichierB_dat.close()
fichierB_dat = open('/net/argos/data/parvati/lahlod/ARCTIC/AMSUB_ice_class/AMSUB89_data_classification_parameters_ice_no-ice_2009.dat','r')
nbtotal=numlines-1
iligne=0
mo = np.zeros([nbtotal],object) # month
spec_emis_thresh = np.zeros([nbtotal],float) # emissivity threshold with spec emis
lamb_emis_thresh = np.zeros([nbtotal],float) # emissivity threshold with spec emis
while (iligne < nbtotal) :
     line = fichierB_dat.readline()
     liste = line.split()
     mo[iligne] = str(liste[0])
     spec_emis_thresh[iligne] = float(liste[7])
     lamb_emis_thresh[iligne] = float(liste[8])
     iligne=iligne+1

fichierB_dat.close()

spec_lim[:] = spec_emis_thresh[vec_read_month] # read spec emis threshold for each month and each frequence
lamb_lim[:] = lamb_emis_thresh[vec_read_month] # read lamb emis threshold for each month and each frequence
for imo in range (0, M):
    print 'month ' + month[imo]
    # AMSUB NETCDF file
    fichierB_nc = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUB_ice_class/cartesian_grid_map_ice_no-ice_' + month[imo] + '2009_AMSUB89_spec_lamb_thresholds.nc', 'r', format='NETCDF3_CLASSIC')
    ice_spec[imo, :, :] = fichierB_nc.variables['spec_ice_area'][:]
    ice_lamb[imo, :, :] = fichierB_nc.variables['lamb_ice_area'][:]
    fichierB_nc.close()
        
    
print 'start mapping'
#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):
    # spec threshold
    print 'map with spec threshold month ' + month[imo]
    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, ice_spec[imo, :, :], -1, 2, 1, colors.ListedColormap(['0.5', 'blue']), 'Sea ice extent with emis spec > ' + str(spec_lim[imo])[0:5])
    title('AMSUB 89GHz - ' + month[imo] + ' 2009')
    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUB/maps/sea_ice_extent/AMSUB89_ice_emis_spec_thresh_' + MONTH[imo] + month[imo] + '2009.png')
    # lamb threshold
    print 'map with lamb threshold month ' + month[imo]
    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, ice_lamb[imo, :, :], -1, 2, 1, colors.ListedColormap(['0.5', 'blue']), 'Sea ice extent with emis lamb > ' + str(lamb_lim[imo])[0:5])
    title('AMSUB 89GHz - ' + month[imo] + ' 2009')
    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUB/maps/sea_ice_extent/AMSUB89_ice_emis_lamb_thresh_' + MONTH[imo] + month[imo] + '2009.png')