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



##########################
# 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 parameters #
###################
x0 = -3000. # min limit of grid
x1 = 2500. # max limit of grid
dx = 10.
xvec = np.arange(x0, x1+dx, dx)
nx = len(xvec) 
y0 = -3000. # min limit of grid
y1 = 3000. # max limit of grid
dy = 10.
yvec = np.arange(y0, y1+dy, dy)
ny = len(yvec)




##########################
# daily OSISAF ice types #
##########################
OSI_TYPE = np.zeros([M, 31, ny, nx], float)
DAY = np.zeros([31], object)
for imo in range (0, M):
    osi_type = np.zeros([month_day[imo], ny, nx], float)
    for ijr in range (0, month_day[imo]):
        if (ijr < 9):
            day = '0' + str(ijr + 1)
        else:
            day = str(ijr + 1)
        DAY[ijr] = day
        print 'date : ' + str(day) + '_' + str(month[imo]) + '_2009'
        osisaf = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/OSISAF_Ice_Types/OSISAF_ice_types_ini/ice_type_nh_polstere-100_multi_2009' + MONTH[imo] + day + '1200.nc', 'r', format='NETCDF3_CLASSIC')
        osi_lat = reshape(osisaf.variables['lat'][:], (851200,1))[:,0]
        osi_lon= reshape(osisaf.variables['lon'][:],(851200,1))[:,0]
        ice_type = reshape(osisaf.variables['ice_type'][:],(851200,1))[:,0]
        osisaf.close()
        bb_osi = nonzero((ice_type > 0))[0]
        print 'start re-gridding' 
        z0 = ice_type[bb_osi].min()
        z1 = ice_type[bb_osi].max()
        zgrid, ngrid, z2grid, sigmagrid, xvec, yvec, xgrid_cart, ygrid_cart = cartesian_grid_monthly.new_cartesian_grid_month(osi_lon[bb_osi], osi_lat[bb_osi], ice_type[bb_osi], z0, z1, dx, dy)
        osi_type[ijr, :, :] = zgrid[:, :]
        osi_lon_grid = xvec
        osi_lat_grid = yvec
    OSI_TYPE[imo, 0 : month_day[imo], :, :] = osi_type[:, :, :]
    print 'stack gridded OSISAF data'
    ################################################
    # stack new-gridded OSISAF data in NETCDF file #
    ################################################
    print 'open file month' + str(month[imo])
    rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/OSISAF_Ice_Types/OSISAF_ice_types_cartesian_grid/OSISAF_ice_types_cartes_grid_res-10_' + 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_lon[:] = xvec
    nc_lat[:] = yvec
    nc_day[:] = np.arange(0., month_day[imo], 1.)
    nc_ice_type[:] = OSI_TYPE[imo, 0 : month_day[imo], :, :]
    rootgrp.close()




########
# test #
########
'''x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50()
x_coast = x_ind
y_coast = y_ind
z_coast = z_ind
zvar = OSI_TYPE[0, 6, :, :]
c0 = 0
c1 = 5
dc = 1
colormap = cm.s3pcpn_l_r
map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, zvar, c0, c1, dc, colormap, 'hello')'''