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


##################################################################################################################
# We devide the loop in two steps : 
# - first loop concerns JANUARY, FEBRUARY, MARCH, APRIL, SEPTEMBER, OCTOBER, NOVEMBER, DECEMBER - use of AMSUA23GHz SPEC emissivity to seperate ice from no-ice zones 
# - second loop concerns MAY, JUNE, JULY, AUGUST - use of AMSUA89GHz SPEC emissivity to seperate ice from no_ice zones
##################################################################################################################


frequ = 89 # apply threshold on this frequency


# daily parameter (2D-array) on ARCTIC area
emis_spec = np.zeros([M, ny, nx, 31], float)
emis_lamb = np.zeros([M, ny, nx, 31], float)
emis_diff = np.zeros([M, ny, nx, 31], float)
emis_ratio = np.zeros([M, ny, nx, 31], float)

# daily parameter (2D-array) on ARCTIC SEA ICE area
daily_spec_ice = np.zeros([M, ny, nx, 31], float)
daily_lamb_ice = np.zeros([M, ny, nx, 31], float)
daily_diff_ice = np.zeros([M, ny, nx, 31], float)
daily_ratio_ice = np.zeros([M, ny, nx, 31], float)


months1 = np.array([0, 1, 2, 3, 8, 9, 10, 11]) # use AMSUA 23GHz to delimit ice/no_ice for JANUARY, FEBRUARY, MARCH, APRIL, SEPTEMBER, OCTOBER, NOVEMBER, DECEMBER
for imo in months1:
    print 'month ' + month[imo]
    ##################################################################################
    # choice of AMSUA 23GHz delimitation ice/no_ice for the sub_classification study #
    ##################################################################################
    print 'open threshold file'
    fichier_emis = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/cartesian_grid_map_ice_no-ice_' + str(month[imo]) + '2009_AMSUA23_spec_lamb_thresholds.nc', 'r', format='NETCDF3_CLASSIC')
    spec_lim = fichier_emis.variables['spec_ice_area'][:] # sea ice pixels defined with spec emis at 23GHz 
    #lamb_lim = fichier_emis.variables['lamb_ice_area'][:]
    fichier_emis.close()
    #########################################################
    # application of AMSUA 23GHz delimitation to other data #
    #########################################################
    print 'open emissivity to sub_classify file'
    ## data of emis SPEC, LAMB, DIFF(SPEC-LAMB)
    fichier_e = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUB' + str(frequ) + '_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC')
    day = fichier_e.variables['days'][:]
    emis_spec[imo, :, :, 0 : month_day[imo]] = fichier_e.variables['e_spec'][:]
    emis_lamb[imo, :, :, 0 : month_day[imo]] = fichier_e.variables['e_lamb'][:]
    emis_diff[imo, :, :, 0 : month_day[imo]] = fichier_e.variables['e_spec_lamb'][:]
    fichier_e.close()
    # calculate emis ratio daily
    for ijr in range (0, month_day[imo]):
        for ilon in range (0, nx):
            for ilat in range (0, ny):
                emis_ratio[imo, ilat, ilon, ijr] = ((emis_lamb[imo, ilat, ilon, ijr] - emis_spec[imo, ilat, ilon, ijr]) / emis_spec[imo, ilat, ilon, ijr]) * 100.
                # create 2D-array of emissivity spec, lamb, diff and ratio on sea ice extent only, defined by AMSUA 23GHz spec emiss threshold   
                if (isnan(spec_lim[ilat, ilon]) == True): # if pixel of sea ice extent defined with spec_emiss_23_threshold corresponds to 'no_ice', then compute NaN in pixel
                    daily_spec_ice[imo, ilat, ilon, ijr] = NaN
                    daily_lamb_ice[imo, ilat, ilon, ijr] = NaN
                    daily_diff_ice[imo, ilat, ilon, ijr] = NaN
                    daily_ratio_ice[imo, ilat, ilon, ijr] = NaN
                else: # if pixel of sea ice extent defined with spec_emiss_23_threshold corresponds to 'ice', then compute value of emis spec, emis lamb or emis diff in pixel
                    daily_spec_ice[imo, ilat, ilon, ijr] = emis_spec[imo, ilat, ilon, ijr]
                    daily_lamb_ice[imo, ilat, ilon, ijr] = emis_lamb[imo, ilat, ilon, ijr]
                    daily_diff_ice[imo, ilat, ilon, ijr] = emis_diff[imo, ilat, ilon, ijr]
                    daily_ratio_ice[imo, ilat, ilon, ijr] = emis_ratio[imo, ilat, ilon, ijr]
    ########################
    # stack in netcdf file #
    ######################## 
    print 'stack in file month ' + str(month[imo])
    rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUB_ice_class/sub_classification/cartesian_grid_map_sea_ice_extent_with-AMSUA23-and-89_' + month[imo] + '2009_AMSUB' + str(frequ) + '_spec_thresholds.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_ice_spec = rootgrp.createVariable('spec_ice_area', 'f', ('latitude', 'longitude', 'days'))
    nc_ice_lamb = rootgrp.createVariable('lamb_ice_area', 'f', ('latitude', 'longitude', 'days'))
    nc_ice_diff = rootgrp.createVariable('diff_ice_area', 'f', ('latitude', 'longitude', 'days'))
    nc_ice_ratio = rootgrp.createVariable('ratio_ice_area', 'f', ('latitude', 'longitude', 'days'))
    nc_lon[:] = xvec
    nc_lat[:] = yvec
    nc_days[:] = np.arange(0, month_day[imo])
    nc_ice_spec[:] = daily_spec_ice[imo, :, :, 0 : month_day[imo]]
    nc_ice_lamb[:] = daily_lamb_ice[imo, :, :, 0 : month_day[imo]]
    nc_ice_diff[:] = daily_diff_ice[imo, :, :, 0 : month_day[imo]]
    nc_ice_ratio[:] = daily_ratio_ice[imo, :, :, 0 : month_day[imo]]
    rootgrp.close()




months2 = np.array([4, 5, 6, 7])# use AMSUA 89GHz to delimit ice/no_ice for MAY, JUNE, JULY, AUGUST
for imo in months2:
    print 'month ' + month[imo]
    ##################################################################################
    # choice of AMSUA 23GHz delimitation ice/no_ice for the sub_classification study #
    ##################################################################################
    print 'open threshold file'
    fichier_emis = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUA_ice_class/cartesian_grid_map_ice_no-ice_' + str(month[imo]) + '2009_AMSUA89_spec_lamb_thresholds.nc', 'r', format='NETCDF3_CLASSIC')
    spec_lim = fichier_emis.variables['spec_ice_area'][:]
    #lamb_lim = fichier_emis.variables['lamb_ice_area'][:]
    fichier_emis.close()
    #########################################################
    # application of AMSUA 23GHz delimitation to other data #
    #########################################################
    print 'open emissivity to sub_classify file'
    ## data of emis SPEC, LAMB, DIFF(SPEC-LAMB)
    fichier_e = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid/res_40/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUB' + str(frequ) + '_' + str(month[imo]) + '2009.nc', 'r', format='NETCDF3_CLASSIC')
    day = fichier_e.variables['days'][:]
    emis_spec[imo, :, :, 0 : month_day[imo]] = fichier_e.variables['e_spec'][:]
    emis_lamb[imo, :, :, 0 : month_day[imo]] = fichier_e.variables['e_lamb'][:]
    emis_diff[imo, :, :, 0 : month_day[imo]] = fichier_e.variables['e_spec_lamb'][:]
    fichier_e.close()
    # calculate emis ratio daily
    for ijr in range (0, month_day[imo]):
        for ilon in range (0, nx):
            for ilat in range (0, ny):
                emis_ratio[imo, ilat, ilon, ijr] = ((emis_lamb[imo, ilat, ilon, ijr] - emis_spec[imo, ilat, ilon, ijr]) / emis_spec[imo, ilat, ilon, ijr]) * 100.
                if (isnan(spec_lim[ilat, ilon]) == True):
                    daily_spec_ice[imo, ilat, ilon, ijr] = NaN
                    daily_lamb_ice[imo, ilat, ilon, ijr] = NaN
                    daily_diff_ice[imo, ilat, ilon, ijr] = NaN
                    daily_ratio_ice[imo, ilat, ilon, ijr] = NaN
                else:
                    daily_spec_ice[imo, ilat, ilon, ijr] = emis_spec[imo, ilat, ilon, ijr]
                    daily_lamb_ice[imo, ilat, ilon, ijr] = emis_lamb[imo, ilat, ilon, ijr]
                    daily_diff_ice[imo, ilat, ilon, ijr] = emis_diff[imo, ilat, ilon, ijr]
                    daily_ratio_ice[imo, ilat, ilon, ijr] = emis_ratio[imo, ilat, ilon, ijr]
    ########################
    # stack in netcdf file #
    ######################## 
    print 'stack in file month ' + str(month[imo])
    rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/AMSUB_ice_class/sub_classification/cartesian_grid_map_sea_ice_extent_with-AMSUA23-and-89_' + month[imo] + '2009_AMSUB' + str(frequ) + '_spec_thresholds.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_ice_spec = rootgrp.createVariable('spec_ice_area', 'f', ('latitude', 'longitude', 'days'))
    nc_ice_lamb = rootgrp.createVariable('lamb_ice_area', 'f', ('latitude', 'longitude', 'days'))
    nc_ice_diff = rootgrp.createVariable('diff_ice_area', 'f', ('latitude', 'longitude', 'days'))
    nc_ice_ratio = rootgrp.createVariable('ratio_ice_area', 'f', ('latitude', 'longitude', 'days'))
    nc_lon[:] = xvec
    nc_lat[:] = yvec
    nc_days[:] = np.arange(0, month_day[imo])
    nc_ice_spec[:] = daily_spec_ice[imo, :, :, 0 : month_day[imo]]
    nc_ice_lamb[:] = daily_lamb_ice[imo, :, :, 0 : month_day[imo]]
    nc_ice_diff[:] = daily_diff_ice[imo, :, :, 0 : month_day[imo]]
    nc_ice_ratio[:] = daily_ratio_ice[imo, :, :, 0 : month_day[imo]]
    rootgrp.close()





'''
# test:
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):
    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, ratio_ice[imo, :, :], 0., 4., 0.01, cm.s3pcpn_l_r, 'Sea ice extent - emissivity RATIO')
    title('AMSUA ' + str(frequ) + ' - ' + str(month[imo]) + ' 2009')
    plt.savefig('/usr/home/lahlod/twice_d/fig_output_ARCTIC/fig_output_sea_ice_study/ice_class_AMSUA/sub_classification/maps_sea_ice_extent/emiss_ratio_map_AMSUA'+str(frequ)+'_with_AMSUA23-and-30-spec_'+str(month[imo])+'_2009.png')
'''