| 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 scipy.special |
|---|
| 11 | import ffgrid2 |
|---|
| 12 | import map_ffgrid |
|---|
| 13 | from matplotlib import colors |
|---|
| 14 | import cartesian_grid_monthly |
|---|
| 15 | import arctic_map |
|---|
| 16 | import map_cartesian_grid |
|---|
| 17 | |
|---|
| 18 | |
|---|
| 19 | |
|---|
| 20 | ########################## |
|---|
| 21 | # time period parameters # |
|---|
| 22 | ########################## |
|---|
| 23 | MONTH = np.array(['01', '02', '03', '04', '05', '06', '07', '08', '09', '10', '11', '12']) |
|---|
| 24 | month = np.array(['JANUARY', 'FEBRUARY', 'MARCH', 'APRIL', 'MAY', 'JUNE', 'JULY', 'AUGUST', 'SEPTEMBER', 'OCTOBER', 'NOVEMBER', 'DECEMBER']) |
|---|
| 25 | month_day = np.array([31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]) |
|---|
| 26 | M = len(month) |
|---|
| 27 | |
|---|
| 28 | |
|---|
| 29 | |
|---|
| 30 | ######################## |
|---|
| 31 | # grid characteristics # |
|---|
| 32 | ######################## |
|---|
| 33 | x0 = -3000. # min limit of grid |
|---|
| 34 | x1 = 2500. # max limit of grid |
|---|
| 35 | dx = 40. |
|---|
| 36 | xvec = np.arange(x0, x1+dx, dx) |
|---|
| 37 | nx = len(xvec) |
|---|
| 38 | y0 = -3000. # min limit of grid |
|---|
| 39 | y1 = 3000. # max limit of grid |
|---|
| 40 | dy = 40. |
|---|
| 41 | yvec = np.arange(y0, y1+dy, dy) |
|---|
| 42 | ny = len(yvec) |
|---|
| 43 | |
|---|
| 44 | |
|---|
| 45 | |
|---|
| 46 | amb_ice = np.zeros([M, ny, nx, 31], float) |
|---|
| 47 | es_a = np.zeros([M, ny, nx, 31], float) |
|---|
| 48 | el_a = np.zeros([M, ny, nx, 31], float) |
|---|
| 49 | esl_a = np.zeros([M, ny, nx, 31], float) |
|---|
| 50 | es_b = np.zeros([M, ny, nx, 31], float) |
|---|
| 51 | el_b = np.zeros([M, ny, nx, 31], float) |
|---|
| 52 | esl_b = np.zeros([M, ny, nx, 31], float) |
|---|
| 53 | for imo in range (0, M): |
|---|
| 54 | ################### |
|---|
| 55 | # Read AMSUA data # |
|---|
| 56 | ################### |
|---|
| 57 | print 'open AMSUA23 file ' + str(month[imo]) |
|---|
| 58 | file_amsua23 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUA23_' + month[imo] + '2009.nc', 'r', format='NETCDF3_CLASSIC') |
|---|
| 59 | lon_amsua = file_amsua23.variables['longitude'][:] |
|---|
| 60 | lat_amsua = file_amsua23.variables['latitude'][:] |
|---|
| 61 | days_amsua = file_amsua23.variables['days'][:] |
|---|
| 62 | es_a[imo, :, :, 0 : month_day[imo]] = file_amsua23.variables['e_spec'][:] |
|---|
| 63 | el_a[imo, :, :, 0 : month_day[imo]] = file_amsua23.variables['e_lamb'][:, :, 0 : month_day[imo]] |
|---|
| 64 | esl_a[imo, :, :, 0 : month_day[imo]] = file_amsua23.variables['e_spec_lamb'][:, :, 0 : month_day[imo]] |
|---|
| 65 | file_amsua23.close() |
|---|
| 66 | ################### |
|---|
| 67 | # Read AMSUB data # |
|---|
| 68 | ################### |
|---|
| 69 | print 'open AMSUB89 file ' + str(month[imo]) |
|---|
| 70 | file_amsub89 = Dataset('/net/argos/data/parvati/lahlod/ARCTIC/monthly_GLACE/gridded_data/cartesian_grid_monthly_data_lamb_spec_near_nadir_AMSUB89_' + month[imo] + '2009.nc', 'r', format='NETCDF3_CLASSIC') |
|---|
| 71 | lon_amsub = file_amsub89.variables['longitude'][:] |
|---|
| 72 | lat_amsub = file_amsub89.variables['latitude'][:] |
|---|
| 73 | days_amsua = file_amsub89.variables['days'][:] |
|---|
| 74 | es_b[imo, :, :, 0 : month_day[imo]] = file_amsub89.variables['e_spec'][:, :, 0 : month_day[imo]] |
|---|
| 75 | el_b[imo, :, :, 0 : month_day[imo]] = file_amsub89.variables['e_lamb'][:, :, 0 : month_day[imo]] |
|---|
| 76 | esl_b[imo, :, :, 0 : month_day[imo]] = file_amsub89.variables['e_spec_lamb'][:, :, 0 : month_day[imo]] |
|---|
| 77 | file_amsub89.close() |
|---|
| 78 | #################### |
|---|
| 79 | # Read OSISAF data # |
|---|
| 80 | #################### |
|---|
| 81 | print 'open OSISAF file ' + str(month[imo]) |
|---|
| 82 | 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') |
|---|
| 83 | lon_osi = file_osisaf.variables['longitude'][:] |
|---|
| 84 | lat_osi = file_osisaf.variables['latitude'][:] |
|---|
| 85 | days_osi = file_osisaf.variables['days'][:] |
|---|
| 86 | ice_type_osi = file_osisaf.variables['osi_ice_type'][:] |
|---|
| 87 | file_osisaf.close() |
|---|
| 88 | ############### |
|---|
| 89 | # daily study # |
|---|
| 90 | ############### |
|---|
| 91 | for ijr in range (0, month_day[imo]): |
|---|
| 92 | print 'day ' + str(ijr) |
|---|
| 93 | for ilat in range (0, len(lat_osi)): |
|---|
| 94 | for ilon in range (0, len(lon_osi)): |
|---|
| 95 | if (ice_type_osi[ijr, ilat, ilon] > 3.): # select only data points located in 'ambiguous ice type' zones defined by OSISAF |
|---|
| 96 | amb_ice[imo, ilat, ilon, ijr] = 1. |
|---|
| 97 | |
|---|
| 98 | |
|---|
| 99 | # select a pixel to study (0, 0, 33, 62) |
|---|
| 100 | |
|---|
| 101 | |
|---|
| 102 | #90, 51 |
|---|
| 103 | #95, 82 |
|---|
| 104 | |
|---|
| 105 | spec_amsua = np.append(es_a[0, 95, 82, 0 : 31], es_a[1, 95, 82, 0 : 28]) |
|---|
| 106 | spec_amsub = np.append(es_b[0, 95, 82, 0 : 31], es_b[1, 95, 82, 0 : 28]) |
|---|
| 107 | ambig = np.append(amb_ice[0, 95, 82, 0 : 31], amb_ice[1, 95, 82, 0 : 28]) |
|---|
| 108 | for imo in range (2, M): |
|---|
| 109 | spec_amsua = np.append(spec_amsua, es_a[imo, 95, 82, 0 : month_day[imo]]) |
|---|
| 110 | spec_amsub = np.append(spec_amsub, es_b[imo, 95, 82, 0 : month_day[imo]]) |
|---|
| 111 | ambig = np.append(ambig, amb_ice[imo, 95, 82, 0 : month_day[imo]]) |
|---|
| 112 | |
|---|
| 113 | |
|---|
| 114 | # seperately AMSUA and AMSUB |
|---|
| 115 | ion() |
|---|
| 116 | figure() |
|---|
| 117 | plot(np.arange(0, 365, 1), spec_amsua, 'b', label = 'AMSUA 23GHz') |
|---|
| 118 | plot(np.arange(0, 365, 1), spec_amsub, 'g', label = 'AMSUB 89GHz') |
|---|
| 119 | plot(np.arange(0, 365, 1)[nonzero(ambig == 1.)], spec_amsua[nonzero(ambig == 1.)], 'ob', label = 'ambiguous area') |
|---|
| 120 | plot(np.arange(0, 365, 1)[nonzero(ambig == 1.)], spec_amsub[nonzero(ambig == 1.)], 'og', label = 'ambiguous area') |
|---|
| 121 | legend(loc = 3) |
|---|
| 122 | xlim(-3, 368) |
|---|
| 123 | xticks(np.array([0, 31, 59, 90, 120, 151, 182, 212, 243, 273, 304, 334]), month, rotation = 20) |
|---|
| 124 | grid() |
|---|
| 125 | |
|---|
| 126 | # bias AMSUA-AMSUB |
|---|
| 127 | figure() |
|---|
| 128 | plot(np.arange(0, 365, 1), spec_amsua - spec_amsub, 'r', label = 'AMSUA23 - AMSUB89') |
|---|
| 129 | plot(np.arange(0, 365, 1)[nonzero(ambig == 1.)], (spec_amsua - spec_amsub)[nonzero(ambig == 1.)], 'or', label = 'ambiguous area') |
|---|
| 130 | plot(np.arange(0, 365, 1), np.zeros([365], float), 'k') |
|---|
| 131 | legend(loc = 3) |
|---|
| 132 | xlim(-3, 368) |
|---|
| 133 | xticks(np.array([0, 31, 59, 90, 120, 151, 182, 212, 243, 273, 304, 335]), month, rotation = 20) |
|---|
| 134 | |
|---|
| 135 | |
|---|
| 136 | # anomaly with a 15 day time range |
|---|
| 137 | spec_anom = np.zeros([365], float) |
|---|
| 138 | spec_diff = spec_amsua - spec_amsub |
|---|
| 139 | for ii in range (15, 365 - 15): |
|---|
| 140 | spec_anom[ii] = (spec_diff[ii]) - mean(spec_diff[ii - 15 : ii + 15][nonzero(isnan(spec_diff[ii - 15 : ii + 15]) == False)]) |
|---|
| 141 | |
|---|
| 142 | figure() |
|---|
| 143 | plot(np.arange(0, 365, 1), spec_anom, 'r') |
|---|
| 144 | plot(np.arange(0, 365, 1)[nonzero(ambig == 1.)], spec_anom[nonzero(ambig == 1.)], 'or', label = 'ambiguous area') |
|---|
| 145 | plot(np.arange(0, 365, 1), np.zeros([365], float), 'k') |
|---|
| 146 | legend(loc = 3) |
|---|
| 147 | xlim(-3, 368) |
|---|
| 148 | xticks(np.array([0, 31, 59, 90, 120, 151, 182, 212, 243, 273, 304, 335]), month, rotation = 20) |
|---|
| 149 | |
|---|
| 150 | |
|---|
| 151 | |
|---|
| 152 | |
|---|
| 153 | |
|---|
| 154 | # plot of anomaly (espec amsub - espec amsua) on a 15 day sliding base (15 day correct to simulate season transition) |
|---|
| 155 | |
|---|
| 156 | colormap = cm.s3pcpn_l_r |
|---|
| 157 | x_ind, y_ind, z_ind, volume = arctic_map.arctic_map_lat50() |
|---|
| 158 | x_coast = x_ind |
|---|
| 159 | y_coast = y_ind |
|---|
| 160 | z_coast = z_ind |
|---|
| 161 | colormap = colors.ListedColormap(['cyan','blue']) |
|---|
| 162 | map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, amb_es_a[9, ijr, :, :], 0.4, 1.02, 0.01, colormap, 'test') |
|---|
| 163 | for imo in range (0, M): |
|---|
| 164 | map_cartesian_grid.draw_map_cartes_l(x_coast, y_coast, z_coast, volume, xvec, yvec, amb_ice[imo, :, :, 5], 0, 1.5, 0.5, colormap, 'ambiguous ice zone') |
|---|
| 165 | title(month[imo] + ' 2009 - OSISAF') |
|---|
| 166 | plt.savefig('/usr/home/lahlod/twice_d/figure_output_ARCTIC/figure_output_CEN/OSISAF_ice_types/cartesian_grid/OSISAF_cartesian_grid_ambiguous_ice_tye_' + month[imo] + '2009.png') |
|---|
| 167 | |
|---|
| 168 | |
|---|
![(please configure the [header_logo] section in trac.ini)](/twice/chrome/site/your_project_logo.png){kind=logo}
