source: trunk/src/scripts_Laura/ARCTIC/Travail_CEN/ice_delimit_OSISAF_dat.py @ 56

Last change on this file since 56 was 54, checked in by lahlod, 10 years ago

modifs
File size: 6.1 KB
Line 
1#!/usr/bin/env python
2# -*- coding: utf-8 -*-
3import string
4import numpy as np
5import matplotlib.pyplot as plt
6from pylab import *
7from mpl_toolkits.basemap import Basemap
8from mpl_toolkits.basemap import shiftgrid, cm
9from netCDF4 import Dataset
10import arctic_map # function to regrid coast limits
11import cartesian_grid_test # function to convert grid from polar to cartesian
12import scipy.special
13import ffgrid2
14import map_ffgrid
15from matplotlib import colors
16from matplotlib.font_manager import FontProperties
17import map_cartesian_grid
18
19MONTH = np.array(['01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12'])
20month = np.array(['JANUARY', 'FEBRUARY', 'MARCH', 'APRIL', 'MAY', 'JUNE', 'JULY', 'AUGUST', 'SEPTEMBER', 'OCTOBER', 'NOVEMBER', 'DECEMBER'])
21month_day = np.array([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31])
22M = len(month)
23
24
25########################
26# grid characteristics #
27########################
28x0 = -3000. # min limit of grid
29x1 = 2500. # max limit of grid
30dx = 40.
31xvec = np.arange(x0, x1+dx, dx)
32nx = len(xvec) 
33y0 = -3000. # min limit of grid
34y1 = 3000. # max limit of grid
35dy = 40.
36yvec = np.arange(y0, y1+dy, dy)
37ny = len(yvec)
38
39
40
41
42#########################################
43# comparison with OSISAF sea ice extent #
44#########################################
45osi_ice = np.zeros([M, 31, 151, 139])
46for imo in range (0, M):
47    print 'open OSISAF file month ' + str(month[imo])
48    fichier_osisaf = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/OSISAF_Ice_Types/OSISAF_ice_types_cartesian_grid/OSISAF_ice_types_cartes_grid_' + str(month[imo]) + '_2009.nc', 'r', format='NETCDF3_CLASSIC')
49    xdist = fichier_osisaf.variables['longitude'][:]
50    ydist = fichier_osisaf.variables['latitude'][:]
51    day = fichier_osisaf.variables['days'][:]
52    osi_type = fichier_osisaf.variables['osi_ice_type'][:]
53    fichier_osisaf.close()
54    for ilon in range (0, nx):
55        for ilat in range (0, ny):
56            for ijr in range (0, month_day[imo]):
57                if ((xdist[ilon] >= -312.5) & (xdist[ilon] <= 328.5) & (ydist[ilat] >=-334.3) & (ydist[ilat] <= 338.5)):
58                    osi_ice[imo, ijr, ilat, ilon] = 1.
59                else:
60                    if (osi_type[ijr, ilat, ilon] > 1.):
61                        osi_ice[imo, ijr, ilat, ilon] = 1.
62                    else:
63                        if (isnan(osi_type[ijr, ilat, ilon]) == True):
64                            osi_ice[imo, ijr, ilat, ilon] = NaN
65                        else: 
66                            osi_ice[imo, ijr, ilat, ilon] = NaN
67
68
69# test:
70ion()
71x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50()
72x_coast = x_ind
73y_coast = y_ind
74z_coast = z_ind
75colormap = colors.ListedColormap(['0.5', 'cyan','blue', 'green'])
76map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, osi_type[1, :, :], 0, 5, 1, cm.s3pcpn_l_r, 'ice_type')
77map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, osi_ice[0, 0, :, :], -1., 2, 1, cm.s3pcpn_l_r, 'test')
78
79
80
81
82##################################
83# calculation of OSISAF ice area #
84##################################
85pix_area = 40. * 40.
86nb_pix_osisaf = np.zeros([M, 31], float)
87total_ice_area_osisaf = np.zeros([M, 31], float)
88total_monthly_ice_area_osisaf = np.zeros([M], float)
89std_total_monthly_ice_area_osisaf = np.zeros([M], float)
90for imo in range (0, M):
91    for ijr in range (0, month_day[imo]):
92        ice = reshape(osi_ice[imo, ijr, :, :], size(osi_ice[imo, ijr, :, :]))
93        nb_pix_osisaf[imo, ijr] = len(np.where(ice == 1.)[0])
94        total_ice_area_osisaf[imo, ijr] = (pix_area * nb_pix_osisaf[imo, ijr]) # daily total ice area
95    total_monthly_ice_area_osisaf[imo] = mean(total_ice_area_osisaf[imo, 0 : month_day[imo]][nonzero(isnan(total_ice_area_osisaf[imo, 0 : month_day[imo]]) == False)]) # monthly mean total ice area
96    alpha = 1. / float(month_day[imo] - 1)
97    area_mean = mean(total_ice_area_osisaf[imo, :][nonzero(isnan(total_ice_area_osisaf[imo, :]) == False)])
98    std_sum = sum((total_ice_area_osisaf[imo, :] - area_mean) ** 2)
99    std_total_monthly_ice_area_osisaf[imo] = sqrt(alpha * std_sum)
100
101
102
103
104
105
106
107
108#################################
109# stack of ice area OSISAF data #
110#################################
111
112#### total ice area in .DAT file ####
113print 'start stacking in .dat file'
114data_osisaf = open ('/net/argos/data/parvati/lahlod/ARCTIC/OSISAF_Ice_Types/OSISAF_daily_ice_no-ice_2009.dat', 'a')
115for imo in range (0, M):
116    for ijr in range (0, month_day[imo]):
117        data_osisaf.write(('%(months)10s    %(days)2.0f    %(total_monthly_ice_area_osisaf)10.5f    %(std_total_monthly_ice_area_osisaf)10.5f    %(total_ice_area_osisaf)10.5f   \n' % {
118                                            'months':month[imo],
119                                            'days':np.arange(1, month_day[imo] + 1)[ijr],
120                                            'total_monthly_ice_area_osisaf':total_monthly_ice_area_osisaf[imo],
121                                            'std_total_monthly_ice_area_osisaf':std_total_monthly_ice_area_osisaf[imo],
122                                            'total_ice_area_osisaf':total_ice_area_osisaf[imo, ijr],
123                                            }))
124
125data_osisaf.close()
126
127
128#### OSISAF sea ice extent maps in NETCDF monthly files ####
129print 'stacking of gridded data'
130for imo in range (0, M):
131    print 'file for month ' + str(month[imo])
132    rootgrp = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/OSISAF_Ice_Types/OSISAF_ice_types_cartesian_grid/OSISAF_daily_ice_no-ice_' + month[imo] + '2009.nc', 'w', format='NETCDF3_CLASSIC')
133    rootgrp.createDimension('longitude', nx)
134    rootgrp.createDimension('latitude', ny)
135    rootgrp.createDimension('days', month_day[imo])
136    nc_x = rootgrp.createVariable('longitude', 'f', ('longitude',))
137    nc_y = rootgrp.createVariable('latitude', 'f', ('latitude',))
138    nc_day = rootgrp.createVariable('days', 'f', ('days',))
139    nc_ice = rootgrp.createVariable('ice_osi', 'f', ('days', 'latitude', 'longitude'))
140    nc_x[:] = xvec
141    nc_y[:] = yvec
142    nc_day[:] = np.arange(0, month_day[imo], 1)
143    nc_ice[:] = osi_ice[imo, 0 : month_day[imo], :, :] # = 0 if open sea // = 1 if sea ice // = NaN if land
144    rootgrp.close()
Note: See TracBrowser for help on using the repository browser.