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



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)



osi_type = np.zeros([M, 31, ny, nx], float) # daily ice type in polar lon/lat grid
spec = np.zeros([M, 31, ny, nx], float) # daily emissivity spec in polar lon/lat grid
lamb = np.zeros([M, 31, ny, nx], float) # daily emissivity lamb in polar lon/lat grid

for imo in range (3, 4):
    print 'month ' + month[imo]
    # osisaf
    print 'read OSISAF'
    fichier_osi = Dataset('/mma/hermozol/Documents/Data/OSISAF_Ice_Types/OSISAF_ice_types_ffgrid/OSISAF_ice_types_ffgrid_' + month[imo] + '_2009.nc', 'r', format='NETCDF3_CLASSIC')
    lon_osi = fichier_osi.variables['longitude'][:]
    lat_osi = fichier_osi.variables['latitude'][:]
    days_osi = fichier_osi.variables['days'][:]
    osi_ice_type = fichier_osi.variables['osi_ice_type'][:]
    osi_type[imo, 0 : month_day[imo], :, :] = osi_ice_type[:, :, :]
    fichier_osi.close()
    # AMSUB data
    print 'read AMSUB'
    fichier_amsub = open('/mma/hermozol/Documents/Data/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009.dat','r')
    numlines = 0
    for line in fichier_amsub: numlines += 1
    fichier_amsub.close
    fichier_amsub = open('/mma/hermozol/Documents/Data/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009.dat','r')  
    nbtotal = numlines-1
    iligne = 0
    jjr = np.zeros([nbtotal],float)
    lat = np.zeros([nbtotal],float)
    lon = np.zeros([nbtotal],float)
    e_spec = np.zeros([nbtotal],float)
    e_lamb = np.zeros([nbtotal],float)
    while (iligne < nbtotal) :
         line=fichier_amsub.readline()
         liste = line.split()
         jjr[iligne] = float(liste[0])
         lat[iligne] = float(liste[2])
         lon[iligne] = float(liste[1])
         e_spec[iligne] = float(liste[16])
         e_lamb[iligne] = float(liste[17])
         iligne=iligne+1
    fichier_amsub.close()
    for ijr in range (0, month_day[imo]): # read daily data
        # stop computing for this month if day > maximum day in month
        if ((ijr + 1) > month_day[imo]): continue
        else:
           bbjr = nonzero(jjr == float(ijr) + 1.)[0]
           if (len(jjr[bbjr]) == 0):
               print 'PAS DE JOUR : ', str(ijr + 1) 
               continue
           else:
               # re-grid in polar lon/lat daily emis spec and emis_lamb 
               print 'date ', jjr[bbjr][0] 
               lat_jr = lat[bbjr]
               lon_jr = lon[bbjr]
               e_spec_jr = e_spec[bbjr]
               e_lamb_jr = e_lamb[bbjr]
               # spec
               print 'grid spec'
               z0 = e_spec_jr[nonzero(isnan(e_spec_jr) == False)].min()
               z1 = e_spec_jr[nonzero(isnan(e_spec_jr) == False)].max()
               zgrid, xvec, yvec = ffgrid2.ffgrid(lon_jr, lat_jr, e_spec_jr, dx, dy, x0, x1, y0, y1, z0, z1)
               spec [imo, ijr, :, :] = zgrid
               # lamb
               print 'grid lamb'
               z0 = e_lamb_jr.min()
               z1 = e_lamb_jr.max()
               zgrid, xvec, yvec = ffgrid2.ffgrid(lon_jr, lat_jr, e_lamb_jr, dx, dy, x0, x1, y0, y1, z0, z1)
               lamb[imo, ijr, :, :] = zgrid



###################################################################
# stacking of daily gridded data AMSUB and OSISAF in NETCDF files #
###################################################################
for imo in range (11, M):
    rootgrp = Dataset('/mma/hermozol/Documents/Data/daily_OSISAF_SPEC_LAMB/daily_OSISAF_emis_spec_lamb_AMSUB_ffgrid_' + str(month[imo]) + '_2009.nc', 'w', format='NETCDF3_CLASSIC')
    rootgrp.createDimension('longitude', len(xvec))
    rootgrp.createDimension('latitude', len(yvec))
    rootgrp.createDimension('days', month_day[imo])
    nc_lon = rootgrp.createVariable('longitude', 'f', ('longitude',))
    nc_lat = rootgrp.createVariable('latitude', 'f', ('latitude',))
    nc_day = rootgrp.createVariable('days', 'f', ('days',))
    nc_ice_type = rootgrp.createVariable('osi_ice_type', 'f', ('days', 'latitude', 'longitude'))
    nc_spec = rootgrp.createVariable('spec_emis', 'f', ('days', 'latitude', 'longitude'))
    nc_lamb = rootgrp.createVariable('lamb_emis', 'f', ('days', 'latitude', 'longitude'))
    nc_lon[:] = xvec
    nc_lat[:] = yvec
    nc_day[:] = np.arange(0., month_day[imo], 1.)
    nc_ice_type[:] = osi_type[imo, 0 : month_day[imo], :, :]
    nc_spec[:] = spec[imo, 0 : month_day[imo], :, :]
    nc_lamb[:] = lamb[imo, 0 : month_day[imo], :, :]
    rootgrp.close()




'''
#########
# tests #
#########
# spec
# lamb
# osi_type
map_ffgrid.draw_map_l (xvec, yvec, 0.4, 1.02, 0.01, spec[11, 12, :, :] , cm.s3pcpn_l_r, 'emis SPEC')
map_ffgrid.draw_map_l (xvec, yvec, 0.4, 1.02, 0.01, lamb[11, 0, :, :] , cm.s3pcpn_l_r, 'emis LAMB')
cmap2 = colors.ListedColormap(['1.', '0.25', '0.50', '0.75'])
map_ffgrid.draw_map_l (xvec, yvec, 0, 5, 1, osi_type[11, 19, :, :] , cmap2, 'OSISAF ice type')
'''