source: trunk/src/scripts_Laura/ARCTIC/Travail_CEN/read_cartesian_osisaf_ice_types.py @ 55

#!/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 scipy.special
import ffgrid2
import map_ffgrid
from matplotlib import colors
import cartesian_grid_monthly
import arctic_map
import map_cartesian_grid




##########################
# time period parameters #
##########################
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 OSISAF data #
####################
ice_type_month = np.zeros([M, ny, nx], float)
for imo in range (0, M):
    print 'open file ' + str(month[imo])
    file_osisaf = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/OSISAF_Ice_Types/OSISAF_ice_types_cartesian_grid/OSISAF_ice_types_cartes_grid_' + month[imo] + '_2009.nc', 'r', format='NETCDF3_CLASSIC')
    lon_osi = file_osisaf.variables['longitude'][:]
    lat_osi = file_osisaf.variables['latitude'][:]
    days = file_osisaf.variables['days'][:]
    ice_type_osi = file_osisaf.variables['osi_ice_type'][:]
    ice_type_month[imo, :, :] = ice_type_osi[0, :, :]
    file_osisaf.close()



###################
# Map OSISAF data #
###################
for imo in range (0, M):
    x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() # definition of the map coast background //se 'arctic_map.py' function for more details
    x_coast = x_ind 
    y_coast = y_ind
    z_coast = z_ind
    colormap = colors.ListedColormap(['0.20', 'cyan','blue','0.75'])
    map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, ice_type_month[imo, :, :], 0, 5, 1, colormap, 'ice_type')
    title(month[imo] + ' 2009 - OSISAF')
    savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/OSISAF_ice_types/cartesian_grid/OSISAF_cartesian_gridded_ice_types_01' + month[imo] + '2009.png')