[40] | 1 | #!/usr/bin/env python |
---|
| 2 | # -*- coding: utf-8 -*- |
---|
| 3 | import string |
---|
| 4 | import numpy as np |
---|
| 5 | import matplotlib.pyplot as plt |
---|
| 6 | from pylab import * |
---|
| 7 | from mpl_toolkits.basemap import Basemap |
---|
| 8 | from mpl_toolkits.basemap import shiftgrid, cm |
---|
| 9 | from netCDF4 import Dataset |
---|
| 10 | import arctic_map # function to regrid coast limits |
---|
| 11 | import cartesian_grid_test # function to convert grid from polar to cartesian |
---|
| 12 | import scipy.special |
---|
| 13 | import ffgrid2 |
---|
| 14 | import map_ffgrid |
---|
| 15 | from matplotlib import colors |
---|
| 16 | |
---|
| 17 | |
---|
| 18 | |
---|
| 19 | MONTH = np.array(['01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12']) |
---|
| 20 | month = np.array(['JANUARY', 'FEBRUARY', 'MARCH', 'APRIL', 'MAY', 'JUNE', 'JULY', 'AUGUST', 'SEPTEMBER', 'OCTOBER', 'NOVEMBER', 'DECEMBER']) |
---|
| 21 | month_day = np.array([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]) |
---|
| 22 | M = len(month) |
---|
| 23 | |
---|
| 24 | |
---|
| 25 | |
---|
| 26 | ################### |
---|
| 27 | # grid parameters # |
---|
| 28 | ################### |
---|
| 29 | dx=1 |
---|
| 30 | dy=0.5 |
---|
| 31 | x0, x1 = -180., 180. |
---|
| 32 | y0, y1 = 30., 90. |
---|
| 33 | xvec = np.arange(x0, x1 + dx, dx) |
---|
| 34 | yvec = np.arange(y0, y1 + dy, dy) |
---|
| 35 | nx = len(xvec) |
---|
| 36 | ny = len(yvec) |
---|
| 37 | |
---|
| 38 | |
---|
| 39 | |
---|
[56] | 40 | osi_type = np.zeros([M, 31, ny, nx], float) # daily ice type in polar lon/lat grid |
---|
| 41 | spec = np.zeros([M, 31, ny, nx], float) # daily emissivity spec in polar lon/lat grid |
---|
| 42 | lamb = np.zeros([M, 31, ny, nx], float) # daily emissivity lamb in polar lon/lat grid |
---|
| 43 | |
---|
[40] | 44 | for imo in range (3, 4): |
---|
| 45 | print 'month ' + month[imo] |
---|
| 46 | # osisaf |
---|
| 47 | print 'read OSISAF' |
---|
| 48 | 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') |
---|
| 49 | lon_osi = fichier_osi.variables['longitude'][:] |
---|
| 50 | lat_osi = fichier_osi.variables['latitude'][:] |
---|
| 51 | days_osi = fichier_osi.variables['days'][:] |
---|
| 52 | osi_ice_type = fichier_osi.variables['osi_ice_type'][:] |
---|
| 53 | osi_type[imo, 0 : month_day[imo], :, :] = osi_ice_type[:, :, :] |
---|
| 54 | fichier_osi.close() |
---|
| 55 | # AMSUB data |
---|
| 56 | print 'read AMSUB' |
---|
| 57 | fichier_amsub = open('/mma/hermozol/Documents/Data/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009.dat','r') |
---|
| 58 | numlines = 0 |
---|
| 59 | for line in fichier_amsub: numlines += 1 |
---|
| 60 | fichier_amsub.close |
---|
| 61 | fichier_amsub = open('/mma/hermozol/Documents/Data/monthly_GLACE/lamb_spec_param_near_nadir_' + month[imo] + '2009.dat','r') |
---|
| 62 | nbtotal = numlines-1 |
---|
| 63 | iligne = 0 |
---|
| 64 | jjr = np.zeros([nbtotal],float) |
---|
| 65 | lat = np.zeros([nbtotal],float) |
---|
| 66 | lon = np.zeros([nbtotal],float) |
---|
| 67 | e_spec = np.zeros([nbtotal],float) |
---|
| 68 | e_lamb = np.zeros([nbtotal],float) |
---|
| 69 | while (iligne < nbtotal) : |
---|
| 70 | line=fichier_amsub.readline() |
---|
| 71 | liste = line.split() |
---|
| 72 | jjr[iligne] = float(liste[0]) |
---|
| 73 | lat[iligne] = float(liste[2]) |
---|
| 74 | lon[iligne] = float(liste[1]) |
---|
| 75 | e_spec[iligne] = float(liste[16]) |
---|
| 76 | e_lamb[iligne] = float(liste[17]) |
---|
| 77 | iligne=iligne+1 |
---|
| 78 | fichier_amsub.close() |
---|
[56] | 79 | for ijr in range (0, month_day[imo]): # read daily data |
---|
| 80 | # stop computing for this month if day > maximum day in month |
---|
[40] | 81 | if ((ijr + 1) > month_day[imo]): continue |
---|
| 82 | else: |
---|
| 83 | bbjr = nonzero(jjr == float(ijr) + 1.)[0] |
---|
| 84 | if (len(jjr[bbjr]) == 0): |
---|
| 85 | print 'PAS DE JOUR : ', str(ijr + 1) |
---|
| 86 | continue |
---|
| 87 | else: |
---|
[56] | 88 | # re-grid in polar lon/lat daily emis spec and emis_lamb |
---|
[40] | 89 | print 'date ', jjr[bbjr][0] |
---|
| 90 | lat_jr = lat[bbjr] |
---|
| 91 | lon_jr = lon[bbjr] |
---|
| 92 | e_spec_jr = e_spec[bbjr] |
---|
| 93 | e_lamb_jr = e_lamb[bbjr] |
---|
| 94 | # spec |
---|
| 95 | print 'grid spec' |
---|
| 96 | z0 = e_spec_jr[nonzero(isnan(e_spec_jr) == False)].min() |
---|
| 97 | z1 = e_spec_jr[nonzero(isnan(e_spec_jr) == False)].max() |
---|
| 98 | zgrid, xvec, yvec = ffgrid2.ffgrid(lon_jr, lat_jr, e_spec_jr, dx, dy, x0, x1, y0, y1, z0, z1) |
---|
| 99 | spec [imo, ijr, :, :] = zgrid |
---|
| 100 | # lamb |
---|
| 101 | print 'grid lamb' |
---|
| 102 | z0 = e_lamb_jr.min() |
---|
| 103 | z1 = e_lamb_jr.max() |
---|
| 104 | zgrid, xvec, yvec = ffgrid2.ffgrid(lon_jr, lat_jr, e_lamb_jr, dx, dy, x0, x1, y0, y1, z0, z1) |
---|
| 105 | lamb[imo, ijr, :, :] = zgrid |
---|
| 106 | |
---|
| 107 | |
---|
| 108 | |
---|
| 109 | ################################################################### |
---|
| 110 | # stacking of daily gridded data AMSUB and OSISAF in NETCDF files # |
---|
| 111 | ################################################################### |
---|
| 112 | for imo in range (11, M): |
---|
| 113 | 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') |
---|
| 114 | rootgrp.createDimension('longitude', len(xvec)) |
---|
| 115 | rootgrp.createDimension('latitude', len(yvec)) |
---|
| 116 | rootgrp.createDimension('days', month_day[imo]) |
---|
| 117 | nc_lon = rootgrp.createVariable('longitude', 'f', ('longitude',)) |
---|
| 118 | nc_lat = rootgrp.createVariable('latitude', 'f', ('latitude',)) |
---|
| 119 | nc_day = rootgrp.createVariable('days', 'f', ('days',)) |
---|
| 120 | nc_ice_type = rootgrp.createVariable('osi_ice_type', 'f', ('days', 'latitude', 'longitude')) |
---|
| 121 | nc_spec = rootgrp.createVariable('spec_emis', 'f', ('days', 'latitude', 'longitude')) |
---|
| 122 | nc_lamb = rootgrp.createVariable('lamb_emis', 'f', ('days', 'latitude', 'longitude')) |
---|
| 123 | nc_lon[:] = xvec |
---|
| 124 | nc_lat[:] = yvec |
---|
| 125 | nc_day[:] = np.arange(0., month_day[imo], 1.) |
---|
| 126 | nc_ice_type[:] = osi_type[imo, 0 : month_day[imo], :, :] |
---|
| 127 | nc_spec[:] = spec[imo, 0 : month_day[imo], :, :] |
---|
| 128 | nc_lamb[:] = lamb[imo, 0 : month_day[imo], :, :] |
---|
| 129 | rootgrp.close() |
---|
| 130 | |
---|
| 131 | |
---|
| 132 | |
---|
| 133 | |
---|
[56] | 134 | ''' |
---|
[40] | 135 | ######### |
---|
| 136 | # tests # |
---|
| 137 | ######### |
---|
| 138 | # spec |
---|
| 139 | # lamb |
---|
| 140 | # osi_type |
---|
| 141 | map_ffgrid.draw_map_l (xvec, yvec, 0.4, 1.02, 0.01, spec[11, 12, :, :] , cm.s3pcpn_l_r, 'emis SPEC') |
---|
| 142 | map_ffgrid.draw_map_l (xvec, yvec, 0.4, 1.02, 0.01, lamb[11, 0, :, :] , cm.s3pcpn_l_r, 'emis LAMB') |
---|
| 143 | cmap2 = colors.ListedColormap(['1.', '0.25', '0.50', '0.75']) |
---|
| 144 | map_ffgrid.draw_map_l (xvec, yvec, 0, 5, 1, osi_type[11, 19, :, :] , cmap2, 'OSISAF ice type') |
---|
[56] | 145 | ''' |
---|
[40] | 146 | |
---|
| 147 | |
---|
| 148 | |
---|
| 149 | |
---|
| 150 | |
---|
| 151 | |
---|
| 152 | |
---|
| 153 | |
---|