From 6c2a7607c6f85bb42e54580763462a145c95a7a8 Mon Sep 17 00:00:00 2001 From: Kaitlin Naughten Date: Wed, 18 Jan 2017 16:45:21 +1100 Subject: [PATCH] Bunch of new things, comments and documentation coming soon --- .gitignore | 1 + adv_amery_tsplots.py | 114 ++-- adv_amery_tsplots_indiv.py | 193 +++++++ adv_frazil.py | 93 ++++ adv_freezingpt_slice.py | 52 +- adv_mld.py | 100 ++++ adv_polynyas.py | 59 ++ adv_sss.py | 94 ++++ adv_thickness.py | 93 ++++ adv_timeseries_volume.py | 4 +- aice_hi_seasonal.py | 174 ++---- cice_vectorplot.py | 46 +- colormaps.py | 1058 ++++++++++++++++++++++++++++++++++++ ice2ocn_fwflux.py | 3 + ice_drift_seasonal.py | 230 ++++++++ nsidc_aice_monthly.py | 39 +- nsidc_aice_seasonal.py | 40 +- sose_seasonal.py | 147 +++++ ssflux_tamura_monthly.py | 209 +++++++ ssflux_tamura_seasonal.py | 233 ++++++++ temp_salt_seasonal.py | 4 +- uv_vectorplot.py | 50 +- zonal_plot.py | 2 +- 23 files changed, 2731 insertions(+), 307 deletions(-) create mode 100644 adv_amery_tsplots_indiv.py create mode 100644 adv_frazil.py create mode 100644 adv_mld.py create mode 100644 adv_polynyas.py create mode 100644 adv_sss.py create mode 100644 adv_thickness.py create mode 100644 colormaps.py create mode 100644 ice_drift_seasonal.py create mode 100644 sose_seasonal.py create mode 100644 ssflux_tamura_monthly.py create mode 100644 ssflux_tamura_seasonal.py diff --git a/.gitignore b/.gitignore index d8976a7..0e8e824 100644 --- a/.gitignore +++ b/.gitignore @@ -2,5 +2,6 @@ *.png *.log *.jnl +*.pdf area_comparison rignot_data \ No newline at end of file diff --git a/adv_amery_tsplots.py b/adv_amery_tsplots.py index a539f93..07d7a20 100644 --- a/adv_amery_tsplots.py +++ b/adv_amery_tsplots.py @@ -3,85 +3,74 @@ from matplotlib.pyplot import * from calc_z import * -# For each advection experiment, plot zonal slices of temperature and salinity -# through 71E (Amery Ice Shelf) at the end of the simulation. def adv_amery_tsplots (): - num_simulations = 6 - # Paths to simulation directories - paths = ['/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/c4_lowdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/c4_highdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/a4_lowdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/a4_highdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/u3_lowdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/u3_highdif/'] - # End of figure names for each simulation - labels = ['_c4_lowdif.png', '_c4_highdif.png', '_a4_lowdif.png', '_a4_highdif.png', '_u3_lowdif.png', '_u3_highdif.png'] - # Name of ocean output file to read - ocn_file = 'ocean_avg_0001.nc' - # Timestep to plot (average over last day in 1992) - tstep = 366 - # Longitude to plot + paths = ['/short/m68/kaa561/advection/u3_lim/', '/short/m68/kaa561/advection/c4_l/'] + labels = [r'a) Temperature ($^{\circ}$C), U3_LIM', r'b) Temperature ($^{\circ}$C), C4_LD', 'c) Salinity (psu), U3_LIM', 'd) Salinity (psu), C4_LD'] + file_tail = 'ocean_avg_31dec.nc' + var_names = ['temp', 'salt'] + tstep = 1 #366 lon0 = 71 - # Deepest depth to plot depth_min = -500 - # Bounds on colour scale for each variable - temp_bounds = [-2, 3] - salt_bounds = [33.8, 34.8] - # Bounds on latitudes to plot + scale_min = [-2, 33.8] + scale_max = [3, 34.8] + scale_ticks = [1, 0.2] lat_min = -72 lat_max = -50 - # Grid parameters theta_s = 4.0 theta_b = 0.9 hc = 40 N = 31 - # Build titles for each variable based on longitude - if lon0 < 0: - temp_title = r'Temperature ($^{\circ}$C) at ' + str(int(round(-lon0))) + r'$^{\circ}$W' - salt_title = r'Salinity (psu) at ' + str(int(round(-lon0))) + r'$^{\circ}$W' - # Edit longitude to be between0 and 360, following ROMS convention - lon0 += 360 - else: - temp_title = r'Temperature ($^{\circ}$C) at ' + str(int(round(lon0))) + r'$^{\circ}$E' - salt_title = r'Salinity (psu) at ' + str(int(round(lon0))) + r'$^{\circ}$E' - - # Loop over simulations - for sim in range(num_simulations): - # Loop over variables - for var_name in ['temp', 'salt']: - # Read variable, sea surface height, and grid variables - id = Dataset(paths[sim] + ocn_file, 'r') - data_3d = id.variables[var_name][tstep-1,:,:-15,:] - zeta = id.variables['zeta'][tstep-1,:-15,:] - if sim == 0 and var_name == 'temp': - # Grid variables are the same for all simulations so we - # only need to read them once - h = id.variables['h'][:-15,:] - zice = id.variables['zice'][:-15,:] - lon_2d = id.variables['lon_rho'][:-15,:] - lat_2d = id.variables['lat_rho'][:-15,:] + fig = figure(figsize=(18,12)) + for sim in range(2): + for var in range(2): + id = Dataset(paths[sim]+file_tail, 'r') + data_3d = id.variables[var_names[var]][tstep-1,:,:,:] + zeta = id.variables['zeta'][tstep-1,:,:] + if sim==0 and var==0: + h = id.variables['h'][:,:] + zice = id.variables['zice'][:,:] + lon_2d = id.variables['lon_rho'][:,:] + lat_2d = id.variables['lat_rho'][:,:] id.close() - # Get a 3D array of z-coordinates z_3d, sc_r, Cs_r = calc_z(h, zice, theta_s, theta_b, hc, N, zeta) - # Interpolate the variable, z, and latitude to lon0 data, z, lat = interp_lon(data_3d, z_3d, lat_2d, lon_2d, lon0) - # Set up colour levels for plotting - if var_name == 'temp': - lev = linspace(temp_bounds[0], temp_bounds[1], num=40) - elif var_name == 'salt': - lev = linspace(salt_bounds[0], salt_bounds[1], num=40) - # Plot - fig = figure(figsize=(12,6)) - contourf(lat, z, data, lev, cmap='jet', extend='both') - colorbar() - if var_name == 'temp': - title(temp_title) - elif var_name == 'salt': - title(salt_title) - xlabel('Latitude') - ylabel('Depth (m)') + ax = fig.add_subplot(2, 2, 2*var+sim+1) + img = pcolor(lat, z, data, vmin=scale_min[var], vmax=scale_max[var], cmap='jet') + title(labels[2*var+sim], fontsize=24) + if var == 1: + xlabel('Latitude', fontsize=16) + if sim == 0: + ylabel('Depth (m)', fontsize=16) xlim([lat_min, lat_max]) ylim([depth_min, 0]) - # Save plot - fig.savefig(var_name + labels[sim]) + if sim == 1: + if var == 0: + cbaxes = fig.add_axes([0.93, 0.575, 0.01, 0.3]) + elif var == 1: + cbaxes = fig.add_axes([0.93, 0.125, 0.01, 0.3]) + cbar = colorbar(img, ticks=arange(scale_min[var], scale_max[var]+scale_ticks[var], scale_ticks[var]), cax=cbaxes, extend='both') + cbar.ax.tick_params(labelsize=14) + lat_ticks = arange(lat_min+2, lat_max+1, 5) + ax.set_xticks(lat_ticks) + lat_labels = [] + for val in lat_ticks: + lat_labels.append(str(int(round(-val))) + r'$^{\circ}$S') + ax.set_xticklabels(lat_labels, fontsize=14) + depth_ticks = range(depth_min, 0+100, 100) + ax.set_yticks(depth_ticks) + depth_labels = [] + for val in depth_ticks: + depth_labels.append(str(int(round(-val)))) + ax.set_yticklabels(depth_labels, fontsize=14) + + suptitle(r'71$^{\circ}$E (Amery Ice Shelf), 31 December', fontsize=30) + #fig.show() + fig.savefig('adv_amery_tsplots.png') + + # Linearly interpolate data, z, and latitude to the specified longitude. @@ -187,7 +176,6 @@ def interp_lon_helper (lon, lon0): return ie, iw, coeff1, coeff2 -# Command-line interface if __name__ == "__main__": adv_amery_tsplots() diff --git a/adv_amery_tsplots_indiv.py b/adv_amery_tsplots_indiv.py new file mode 100644 index 0000000..e3955f8 --- /dev/null +++ b/adv_amery_tsplots_indiv.py @@ -0,0 +1,193 @@ +from netCDF4 import Dataset +from numpy import * +from matplotlib.pyplot import * +from calc_z import * + +# For each advection experiment, plot zonal slices of temperature and salinity +# through 71E (Amery Ice Shelf) at the end of the simulation. +def adv_amery_tsplots_indiv (): + + num_simulations = 6 + # Paths to simulation directories + paths = ['/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/c4_lowdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/c4_highdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/a4_lowdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/a4_highdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/u3_lowdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/u3limiters_lowdif/'] + # End of figure names for each simulation + labels = ['_c4_lowdif.png', '_c4_highdif.png', '_a4_lowdif.png', '_a4_highdif.png', '_u3.png', '_u3_lim.png'] + # Name of ocean output file to read + ocn_file = 'ocean_avg_0001.nc' + # Timestep to plot (average over last day in 1992) + tstep = 366 + # Longitude to plot + lon0 = 71 + # Deepest depth to plot + depth_min = -500 + # Bounds on colour scale for each variable + temp_bounds = [-2, 3] + salt_bounds = [33.8, 34.8] + # Bounds on latitudes to plot + lat_min = -72 + lat_max = -50 + + # Grid parameters + theta_s = 4.0 + theta_b = 0.9 + hc = 40 + N = 31 + + # Build titles for each variable based on longitude + if lon0 < 0: + temp_title = r'Temperature ($^{\circ}$C) at ' + str(int(round(-lon0))) + r'$^{\circ}$W' + salt_title = r'Salinity (psu) at ' + str(int(round(-lon0))) + r'$^{\circ}$W' + # Edit longitude to be between0 and 360, following ROMS convention + lon0 += 360 + else: + temp_title = r'Temperature ($^{\circ}$C) at ' + str(int(round(lon0))) + r'$^{\circ}$E' + salt_title = r'Salinity (psu) at ' + str(int(round(lon0))) + r'$^{\circ}$E' + + # Loop over simulations + for sim in range(num_simulations): + # Loop over variables + for var_name in ['temp', 'salt']: + # Read variable, sea surface height, and grid variables + id = Dataset(paths[sim] + ocn_file, 'r') + data_3d = id.variables[var_name][tstep-1,:,:-15,:] + zeta = id.variables['zeta'][tstep-1,:-15,:] + if sim == 0 and var_name == 'temp': + # Grid variables are the same for all simulations so we + # only need to read them once + h = id.variables['h'][:-15,:] + zice = id.variables['zice'][:-15,:] + lon_2d = id.variables['lon_rho'][:-15,:] + lat_2d = id.variables['lat_rho'][:-15,:] + id.close() + # Get a 3D array of z-coordinates + z_3d, sc_r, Cs_r = calc_z(h, zice, theta_s, theta_b, hc, N, zeta) + # Interpolate the variable, z, and latitude to lon0 + data, z, lat = interp_lon(data_3d, z_3d, lat_2d, lon_2d, lon0) + # Set up colour levels for plotting + if var_name == 'temp': + lev = linspace(temp_bounds[0], temp_bounds[1], num=40) + elif var_name == 'salt': + lev = linspace(salt_bounds[0], salt_bounds[1], num=40) + # Plot + fig = figure(figsize=(12,6)) + contourf(lat, z, data, lev, cmap='jet', extend='both') + colorbar() + if var_name == 'temp': + title(temp_title) + elif var_name == 'salt': + title(salt_title) + xlabel('Latitude') + ylabel('Depth (m)') + xlim([lat_min, lat_max]) + ylim([depth_min, 0]) + # Save plot + fig.savefig(var_name + labels[sim]) + + +# Linearly interpolate data, z, and latitude to the specified longitude. +# Input: +# data_3d = array of data, dimension depth x lat x lon +# z_3d = array of depth values (negative, in metres), dimension depth x lat x lon +# lat_2d = array of latitudevalues, dimension lat x lon +# lon_2d = array of longitude values, dimension lat x lon (between -180 and 180) +# lon0 = longitude to interpolate to (between -180 and 180) +# Output: +# data = array of data interpolated to lon0, dimension depth x lat +# z = array of depth values interpolated to lon0, dimension depth x lat +# lat = array of latitude values interpolated to lon0, dimension depth x lat +def interp_lon (data_3d, z_3d, lat_2d, lon_2d, lon0): + + # Save dimensions + num_depth = size(data_3d, 0) + num_lat = size(data_3d, 1) + num_lon = size(data_3d, 2) + # Set up output arrays + data = ma.empty([num_depth, num_lat]) + z = ma.empty([num_depth, num_lat]) + lat = ma.empty([num_depth, num_lat]) + + # Loop over latitudes; can't find a cleaner way to do this + for j in range(num_lat): + # Extract the longitude values of this slice + lon_tmp = lon_2d[j,:] + # Get indices and coefficients for interpolation + ie, iw, coeffe, coeffw = interp_lon_helper(lon_tmp, lon0) + data[:,j] = coeffe*data_3d[:,j,ie] + coeffw*data_3d[:,j,iw] + z[:,j] = coeffe*z_3d[:,j,ie] + coeffw*z_3d[:,j,iw] + lat[:,j] = coeffe*lat_2d[j,ie] + coeffw*lat_2d[j,iw] + + return data, z, lat + + +# Calculate indices and coefficients for linear interpolation of longitude. +# This takes care of all the mod 360 nonsense. +# Input: +# lon = 1D array of longitude values (straight out of ROMS i.e. between slightly < 0 and slightly > 360) +# lon0 = longitude to interpolate to (between 0 and 360) +# Output: +# ie, iw, coeffe, coeffw = integers (ie and iw) and coefficients (coeffe and coeffw) such that +# coeffe*lon[ie] + coeffw*lon[iw] = lon0, which will also hold for any +# variable on this longitude grid. ie is the index of the nearest point +# to the east of lon0; iw the nearest point to the west. +def interp_lon_helper (lon, lon0): + + if lon0 < amin(lon) or lon0 > amax(lon): + # Special case: lon0 on periodic boundary + # Be careful with mod 360 here + + # Find the periodic boundary + dlon = lon[1:] - lon[0:-1] + bdry = argmax(abs(dlon)) + if dlon[bdry] < -300: + # Jumps from almost 360 to just over 0 + iw = bdry + ie = bdry + 1 + else: + # Periodic boundary lines up with the array boundary + iw = size(lon) - 1 + ie = 0 + # Calculate difference between lon0 and lon[iw], mod 360 if necessary + dlon_num = lon0 - lon[iw] + if dlon_num < -300: + dlon_num += 360 + # Calculate difference between lon[ie] and lon[iw], mod 360 + dlon_den = lon[ie] - lon[iw] + 360 + + else: + # General case + + # Add or subtract 360 from longitude values which wrap around + # so that longitude increases monotonically from west to east + i = arange(1, size(lon)+1) + index1 = nonzero((i > 1200)*(lon < 100)) + lon[index1] = lon[index1] + 360 + index2 = nonzero((i < 200)*(lon > 300)) + lon[index2] = lon[index2] - 360 + + # Take mod 360 of lon0 if necessary + if all(lon < lon0): + lon0 -= 360 + if all(lon > lon0): + lon0 += 360 + + # Find the first index eastward of lon0 + ie = nonzero(lon > lon0)[0][0] + # The index before it will be the last index westward of lon0 + iw = ie - 1 + + dlon_num = lon0 - lon[iw] + dlon_den = lon[ie] - lon[iw] + + if dlon_num > 5 or dlon_den > 5: + print 'interp_lon_helper: Problem at periodic boundary' + return + coeff1 = dlon_num/dlon_den + coeff2 = 1 - coeff1 + + return ie, iw, coeff1, coeff2 + + +# Command-line interface +if __name__ == "__main__": + + adv_amery_tsplots_indiv() diff --git a/adv_frazil.py b/adv_frazil.py new file mode 100644 index 0000000..6dca7a2 --- /dev/null +++ b/adv_frazil.py @@ -0,0 +1,93 @@ +from netCDF4 import * +from numpy import * +from matplotlib.pyplot import * +#import colormaps as cmaps + +def adv_frazil (): + + paths = ['/short/m68/kaa561/advection/u3_lim/', '/short/m68/kaa561/advection/u3/', '/short/m68/kaa561/advection/c4_l/', '/short/m68/kaa561/advection/c4_h/', '/short/m68/kaa561/advection/a4_l/', '/short/m68/kaa561/advection/a4_h/'] + labels = ['a) U3_LIM', 'b) U3 - U3_LIM', 'c) C4_LD - U3_LIM', 'd) C4_HD - U3_LIM', 'e) A4_LD - U3_LIM', 'f) A4_HD - U3_LIM'] + file_tail = 'iceh_avg.nc' + + max_abs = 0.2 + tick_abs = 0.05 + max_anom = 1.0 + tick_anom = 0.5 + + deg2rad = pi/180. + lon_c = 50 + lat_c = -83 + radius = 10.5 + circle_bdry = -70+90 + + id = Dataset(paths[0]+file_tail, 'r') + data_tmp = id.variables['frazil'][0,:350,:] + lon_tmp = id.variables['TLON'][:350,:] + lat_tmp = id.variables['TLAT'][:350,:] + mask_tmp = id.variables['tmask'][:350,:] + id.close() + + lon = ma.empty([size(lon_tmp,0), size(lon_tmp,1)+1]) + lat = ma.empty([size(lat_tmp,0), size(lat_tmp,1)+1]) + mask = ma.empty([size(mask_tmp,0), size(mask_tmp,1)+1]) + data0 = ma.empty([size(data_tmp,0), size(data_tmp,1)+1]) + lon[:,:-1] = lon_tmp + lon[:,-1] = lon_tmp[:,0] + lat[:,:-1] = lat_tmp + lat[:,-1] = lat_tmp[:,0] + mask[:,:-1] = mask_tmp + mask[:,-1] = mask_tmp[:,0] + data0[:,:-1] = data_tmp + data0[:,-1] = data_tmp[:,0] + + land = ma.masked_where(mask==1, mask) + + x = -(lat+90)*cos(lon*deg2rad+pi/2) + y = (lat+90)*sin(lon*deg2rad+pi/2) + x_c = -(lat_c+90)*cos(lon_c*deg2rad+pi/2) + y_c = (lat_c+90)*sin(lon_c*deg2rad+pi/2) + x_reg, y_reg = meshgrid(linspace(-circle_bdry, circle_bdry, num=100), linspace(-circle_bdry, circle_bdry, num=100)) + land_circle = zeros(shape(x_reg)) + land_circle = ma.masked_where(sqrt((x_reg-x_c)**2 + (y_reg-y_c)**2) > radius, land_circle) + + fig = figure(figsize=(9,15)) + ax = fig.add_subplot(3, 2, 1, aspect='equal') + contourf(x, y, land, 1, colors=(('0.6', '0.6', '0.6'))) + contourf(x_reg, y_reg, land_circle, 1, colors=(('0.6', '0.6', '0.6'))) + img0 = pcolor(x, y, data0, vmin=0, vmax=max_abs, cmap='jet') + axis('off') + title(labels[0], fontsize=20) + cbaxes0 = fig.add_axes([0.025, 0.7, 0.02, 0.2]) + cbar0 = colorbar(img0, ticks=arange(0, max_abs+tick_abs, tick_abs), cax=cbaxes0, extend='max') + cbar0.ax.tick_params(labelsize=16) + + for sim in range(1, len(paths)): + id = Dataset(paths[sim]+file_tail, 'r') + data_tmp = id.variables['frazil'][0,:350,:] + id.close() + data = ma.empty([size(data_tmp,0), size(data_tmp,1)+1]) + data[:,:-1] = data_tmp + data[:,-1] = data_tmp[:,0] + data = data - data0 + ax = fig.add_subplot(3, 2, sim+1, aspect='equal') + contourf(x, y, land, 1, colors=(('0.6', '0.6', '0.6'))) + contourf(x_reg, y_reg, land_circle, 1, colors=(('0.6', '0.6', '0.6'))) + img = pcolor(x, y, data, vmin=-max_anom, vmax=max_anom, cmap='RdBu_r') + axis('off') + title(labels[sim], fontsize=20) + if sim == 3: + cbaxes = fig.add_axes([0.025, 0.4, 0.02, 0.2]) + cbar = colorbar(img, ticks=arange(-max_anom, max_anom+tick_anom, tick_anom), cax=cbaxes, extend='both') + cbar.ax.tick_params(labelsize=16) + + suptitle('Annually averaged frazil ice formation (cm/day)', fontsize=24) + subplots_adjust(wspace=0.025,hspace=0.2) + + #fig.show() + fig.savefig('adv_frazil.png') + + +if __name__ == '__main__': + + adv_frazil() + diff --git a/adv_freezingpt_slice.py b/adv_freezingpt_slice.py index 309ab26..1c6ca85 100644 --- a/adv_freezingpt_slice.py +++ b/adv_freezingpt_slice.py @@ -10,15 +10,16 @@ def adv_freezingpt_slice (): # Path to ocean history file - file_path = '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/c4_lowdif/ocean_his_0001.nc' + file_path = '/short/m68/kaa561/advection/c4_l/ocean_his_6july.nc' # Timestep to plot - tstep = 189 + tstep = 1 #188 # i-index to plot (1-based) i_val = 1250 # Deepest depth to plot depth_min = -100 # Bounds on colour scale - colour_bounds = [-0.3, 0.3] + scale_max = 0.3 + scale_tick = 0.1 # Bounds on latitudes to plot lat_min = -78 lat_max = -72 @@ -44,7 +45,7 @@ def adv_freezingpt_slice (): id.close() # Calculate freezing point as seen by supercooling code - tfr = -0.054*salt + tfr = salt/(-18.48 + salt*18.48/1000.0) # Calculate difference from freezing point deltat = temp - tfr @@ -54,33 +55,30 @@ def adv_freezingpt_slice (): z = z_3d[:,:,i_val-1] lat = tile(lat_2d[:,i_val-1], (N,1)) - # Determine colour bounds - if colour_bounds is not None: - # Specified by user - scale_min = colour_bounds[0] - scale_max = colour_bounds[1] - if scale_min == -scale_max: - # Centered on zero; use a red-yellow-blue colour scale - colour_map = 'RdYlBu_r' - else: - # Use a rainbow colour scale - colour_map = 'jet' - else: - # Determine automatically - scale_min = amin(deltat) - scale_max = amax(deltat) - colour_map = 'jet' - # Plot (pcolor not contour to show what each individual cell is doing) - fig = figure(figsize=(12,6)) - pcolor(lat, z, deltat, vmin=scale_min, vmax=scale_max, cmap=colour_map) - colorbar() - title(r'Difference from freezing point ($^{\circ}$C) in Weddell Sea, 7 July') - xlabel('Latitude') - ylabel('Depth (m)') + fig, ax = subplots(figsize=(12,6)) + pcolor(lat, z, deltat, vmin=-scale_max, vmax=scale_max, cmap='RdBu_r') + cbar = colorbar(extend='both', ticks=arange(-scale_max, scale_max+scale_tick, scale_tick)) + cbar.ax.tick_params(labelsize=14) + title(r'Difference from freezing point ($^{\circ}$C) in Weddell Sea: C4_LD', fontsize=18) + xlabel('Latitude', fontsize=16) + ylabel('Depth (m)', fontsize=16) xlim([lat_min, lat_max]) ylim([depth_min, 0]) + lat_ticks = arange(lat_min+1, lat_max+1, 1) + xticks(lat_ticks) + lat_labels = [] + for val in lat_ticks: + lat_labels.append(str(int(round(-val))) + r'$^{\circ}$S') + ax.set_xticklabels(lat_labels, fontsize=14) + depth_ticks = arange(depth_min, 0+25, 25) + yticks(depth_ticks) + depth_labels = [] + for val in depth_ticks: + depth_labels.append(str(int(round(-val)))) + ax.set_yticklabels(depth_labels, fontsize=14) + # Finished if save: fig.savefig(fig_name) diff --git a/adv_mld.py b/adv_mld.py new file mode 100644 index 0000000..2fef7bd --- /dev/null +++ b/adv_mld.py @@ -0,0 +1,100 @@ +from netCDF4 import * +from numpy import * +from matplotlib.pyplot import * +#import colormaps as cmaps + +def adv_mld (): + + paths = ['/short/m68/kaa561/advection/u3_lim/', '/short/m68/kaa561/advection/u3/', '/short/m68/kaa561/advection/c4_l/', '/short/m68/kaa561/advection/c4_h/', '/short/m68/kaa561/advection/a4_l/', '/short/m68/kaa561/advection/a4_h/'] + labels = ['a) U3_LIM', 'b) U3 - U3_LIM', 'c) C4_LD - U3_LIM', 'd) C4_HD - U3_LIM', 'e) A4_LD - U3_LIM', 'f) A4_HD - U3_LIM'] + file_tail = 'ocean_avg_23aug.nc' + tstep = 1 #236 + + max_abs = 300 + tick_abs = 100 + max_anom = 200 + tick_anom = 100 + + deg2rad = pi/180. + lon_c = 50 + lat_c = -83 + radius = 10.5 + circle_bdry = -70+90 + + id = Dataset(paths[0]+file_tail, 'r') + data_tmp = id.variables['Hsbl'][tstep-1,:350,:] + lon_tmp = id.variables['lon_rho'][:350,:] + lat_tmp = id.variables['lat_rho'][:350,:] + mask_tmp = id.variables['mask_rho'][:350,:] + zice_tmp = id.variables['zice'][:350,:] + id.close() + + index = zice_tmp != 0 + mask_tmp[index] = 0.0 + data_tmp = ma.masked_where(zice_tmp!=0, -data_tmp) + + lon = ma.empty([size(lon_tmp,0), size(lon_tmp,1)+1]) + lat = ma.empty([size(lat_tmp,0), size(lat_tmp,1)+1]) + mask = ma.empty([size(mask_tmp,0), size(mask_tmp,1)+1]) + data0 = ma.empty([size(data_tmp,0), size(data_tmp,1)+1]) + lon[:,:-1] = lon_tmp + lon[:,-1] = lon_tmp[:,0] + lat[:,:-1] = lat_tmp + lat[:,-1] = lat_tmp[:,0] + mask[:,:-1] = mask_tmp + mask[:,-1] = mask_tmp[:,0] + data0[:,:-1] = data_tmp + data0[:,-1] = data_tmp[:,0] + + land = ma.masked_where(mask==1, mask) + + x = -(lat+90)*cos(lon*deg2rad+pi/2) + y = (lat+90)*sin(lon*deg2rad+pi/2) + x_c = -(lat_c+90)*cos(lon_c*deg2rad+pi/2) + y_c = (lat_c+90)*sin(lon_c*deg2rad+pi/2) + x_reg, y_reg = meshgrid(linspace(-circle_bdry, circle_bdry, num=100), linspace(-circle_bdry, circle_bdry, num=100)) + land_circle = zeros(shape(x_reg)) + land_circle = ma.masked_where(sqrt((x_reg-x_c)**2 + (y_reg-y_c)**2) > radius, land_circle) + + fig = figure(figsize=(9,15)) + ax = fig.add_subplot(3, 2, 1, aspect='equal') + contourf(x, y, land, 1, colors=(('0.6', '0.6', '0.6'))) + contourf(x_reg, y_reg, land_circle, 1, colors=(('0.6', '0.6', '0.6'))) + img0 = pcolor(x, y, data0, vmin=0, vmax=max_abs, cmap='jet') #cmaps.viridis) + axis('off') + title(labels[0], fontsize=20) + cbaxes0 = fig.add_axes([0.025, 0.7, 0.02, 0.2]) + cbar0 = colorbar(img0, ticks=arange(0, max_abs+tick_abs, tick_abs), cax=cbaxes0, extend='max') + cbar0.ax.tick_params(labelsize=16) + + for sim in range(1, len(paths)): + id = Dataset(paths[sim]+file_tail, 'r') + data_tmp = id.variables['Hsbl'][tstep-1,:350,:] + id.close() + data_tmp = ma.masked_where(zice_tmp!=0, -data_tmp) + data = ma.empty([size(data_tmp,0), size(data_tmp,1)+1]) + data[:,:-1] = data_tmp + data[:,-1] = data_tmp[:,0] + data = data - data0 + ax = fig.add_subplot(3, 2, sim+1, aspect='equal') + contourf(x, y, land, 1, colors=(('0.6', '0.6', '0.6'))) + contourf(x_reg, y_reg, land_circle, 1, colors=(('0.6', '0.6', '0.6'))) + img = pcolor(x, y, data, vmin=-max_anom, vmax=max_anom, cmap='RdBu_r') + axis('off') + title(labels[sim], fontsize=20) + if sim == 3: + cbaxes = fig.add_axes([0.025, 0.4, 0.02, 0.2]) + cbar = colorbar(img, ticks=arange(-max_anom, max_anom+tick_anom, tick_anom), cax=cbaxes, extend='both') + cbar.ax.tick_params(labelsize=16) + + suptitle('Mixed layer depth (m) on 23 August', fontsize=28) + subplots_adjust(wspace=0.025,hspace=0.15) + + #fig.show() + fig.savefig('adv_mld.png') + + +if __name__ == '__main__': + + adv_mld() + diff --git a/adv_polynyas.py b/adv_polynyas.py new file mode 100644 index 0000000..956d74a --- /dev/null +++ b/adv_polynyas.py @@ -0,0 +1,59 @@ +from netCDF4 import * +from numpy import * +from matplotlib.pyplot import * + +def adv_polynyas (): + + paths = ['/short/m68/kaa561/advection/u3_lim/', '/short/m68/kaa561/advection/c4_l/'] + labels = ['a) U3_LIM', 'b) C4_LD'] + file_tail = 'iceh.1992-08-23.nc' + lon_min = 100 + lon_max = 140 + lat_min = -67.1 + lat_max = -64.9 + lat_min_label = -67 + lat_max_label = -65 + + fig = figure(figsize=(18,6)) + for sim in range(2): + id = Dataset(paths[sim]+file_tail, 'r') + data = id.variables['aice'][0,70:280,350:750] + if sim == 0: + lon = id.variables['TLON'][70:280,350:750] + lat = id.variables['TLAT'][70:280,350:750] + id.close() + ax = fig.add_subplot(1, 2, sim+1) + img = pcolor(lon, lat, data, vmin=0, vmax=1, cmap='jet') + title(labels[sim], fontsize=20) + xlabel('Longitude', fontsize=16) + ylabel('Latitude', fontsize=16) + xlim([lon_min, lon_max]) + ylim([lat_min, lat_max]) + if sim == 1: + cbaxes = fig.add_axes([0.93, 0.25, 0.015, 0.5]) + cbar = colorbar(ticks=arange(0, 1+0.25, 0.25), cax=cbaxes) + cbar.ax.tick_params(labelsize=14) + + lon_ticks = arange(lon_min, lon_max+10, 10) + ax.set_xticks(lon_ticks) + lon_labels = [] + for val in lon_ticks: + lon_labels.append(str(int(round(val))) + r'$^{\circ}$E') + ax.set_xticklabels(lon_labels, fontsize=14) + lat_ticks = arange(lat_min_label, lat_max_label+1, 1) + ax.set_yticks(lat_ticks) + lat_labels = [] + for val in lat_ticks: + lat_labels.append(str(int(round(-val))) + r'$^{\circ}$S') + ax.set_yticklabels(lat_labels, fontsize=14) + suptitle('Sea ice concentration on 23 August', fontsize=22) + + #fig.show() + fig.savefig('adv_polynyas.png') + + +if __name__ == "__main__": + + adv_polynyas() + + diff --git a/adv_sss.py b/adv_sss.py new file mode 100644 index 0000000..6333b25 --- /dev/null +++ b/adv_sss.py @@ -0,0 +1,94 @@ +from netCDF4 import * +from numpy import * +from matplotlib.pyplot import * +#import colormaps as cmaps + +def adv_sss (): + + paths = ['/short/m68/kaa561/advection/u3_lim/', '/short/m68/kaa561/advection/u3/', '/short/m68/kaa561/advection/c4_l/', '/short/m68/kaa561/advection/c4_h/', '/short/m68/kaa561/advection/a4_l/', '/short/m68/kaa561/advection/a4_h/'] + labels = ['a) U3_LIM', 'b) U3 - U3_LIM', 'c) C4_LD - U3_LIM', 'd) C4_HD - U3_LIM', 'e) A4_LD - U3_LIM', 'f) A4_HD - U3_LIM'] + file_tail = 'iceh.1992-08-23.nc' + + min_abs = 32.75 + max_abs = 34.75 + tick_abs = 0.5 + max_anom = 3.0 + tick_anom = 1.5 + + deg2rad = pi/180. + lon_c = 50 + lat_c = -83 + radius = 10.5 + circle_bdry = -70+90 + + id = Dataset(paths[0]+file_tail, 'r') + data_tmp = id.variables['sss'][0,:350,:] + lon_tmp = id.variables['TLON'][:350,:] + lat_tmp = id.variables['TLAT'][:350,:] + mask_tmp = id.variables['tmask'][:350,:] + id.close() + + lon = ma.empty([size(lon_tmp,0), size(lon_tmp,1)+1]) + lat = ma.empty([size(lat_tmp,0), size(lat_tmp,1)+1]) + mask = ma.empty([size(mask_tmp,0), size(mask_tmp,1)+1]) + data0 = ma.empty([size(data_tmp,0), size(data_tmp,1)+1]) + lon[:,:-1] = lon_tmp + lon[:,-1] = lon_tmp[:,0] + lat[:,:-1] = lat_tmp + lat[:,-1] = lat_tmp[:,0] + mask[:,:-1] = mask_tmp + mask[:,-1] = mask_tmp[:,0] + data0[:,:-1] = data_tmp + data0[:,-1] = data_tmp[:,0] + + land = ma.masked_where(mask==1, mask) + + x = -(lat+90)*cos(lon*deg2rad+pi/2) + y = (lat+90)*sin(lon*deg2rad+pi/2) + x_c = -(lat_c+90)*cos(lon_c*deg2rad+pi/2) + y_c = (lat_c+90)*sin(lon_c*deg2rad+pi/2) + x_reg, y_reg = meshgrid(linspace(-circle_bdry, circle_bdry, num=100), linspace(-circle_bdry, circle_bdry, num=100)) + land_circle = zeros(shape(x_reg)) + land_circle = ma.masked_where(sqrt((x_reg-x_c)**2 + (y_reg-y_c)**2) > radius, land_circle) + + fig = figure(figsize=(9,15)) + ax = fig.add_subplot(3, 2, 1, aspect='equal') + contourf(x, y, land, 1, colors=(('0.6', '0.6', '0.6'))) + contourf(x_reg, y_reg, land_circle, 1, colors=(('0.6', '0.6', '0.6'))) + img0 = pcolor(x, y, data0, vmin=min_abs, vmax=max_abs, cmap='jet') #cmaps.viridis) + axis('off') + title(labels[0], fontsize=20) + cbaxes0 = fig.add_axes([0.025, 0.7, 0.02, 0.2]) + cbar0 = colorbar(img0, ticks=arange(min_abs, max_abs+tick_abs, tick_abs), cax=cbaxes0, extend='both') + cbar0.ax.tick_params(labelsize=16) + + for sim in range(1, len(paths)): + id = Dataset(paths[sim]+file_tail, 'r') + data_tmp = id.variables['sss'][0,:350,:] + id.close() + data = ma.empty([size(data_tmp,0), size(data_tmp,1)+1]) + data[:,:-1] = data_tmp + data[:,-1] = data_tmp[:,0] + data = data - data0 + ax = fig.add_subplot(3, 2, sim+1, aspect='equal') + contourf(x, y, land, 1, colors=(('0.6', '0.6', '0.6'))) + contourf(x_reg, y_reg, land_circle, 1, colors=(('0.6', '0.6', '0.6'))) + img = pcolor(x, y, data, vmin=-max_anom, vmax=max_anom, cmap='RdBu_r') + axis('off') + title(labels[sim], fontsize=20) + if sim == 3: + cbaxes = fig.add_axes([0.025, 0.4, 0.02, 0.2]) + cbar = colorbar(img, ticks=arange(-max_anom, max_anom+tick_anom, tick_anom), cax=cbaxes, extend='both') + cbar.ax.tick_params(labelsize=16) + + suptitle('Sea surface salinity (psu) on 23 August', fontsize=28) + subplots_adjust(wspace=0.025,hspace=0.15) + + #fig.show() + fig.savefig('adv_sss.png') + + +if __name__ == '__main__': + + adv_sss() + diff --git a/adv_thickness.py b/adv_thickness.py new file mode 100644 index 0000000..f9b468f --- /dev/null +++ b/adv_thickness.py @@ -0,0 +1,93 @@ +from netCDF4 import * +from numpy import * +from matplotlib.pyplot import * +#import colormaps as cmaps + +def adv_thickness (): + + paths = ['/short/m68/kaa561/advection/u3_lim/', '/short/m68/kaa561/advection/u3/', '/short/m68/kaa561/advection/c4_l/', '/short/m68/kaa561/advection/c4_h/', '/short/m68/kaa561/advection/a4_l/', '/short/m68/kaa561/advection/a4_h/'] + labels = ['a) U3_LIM', 'b) U3 - U3_LIM', 'c) C4_LD - U3_LIM', 'd) C4_HD - U3_LIM', 'e) A4_LD - U3_LIM', 'f) A4_HD - U3_LIM'] + file_tail = 'iceh.1992-08-23.nc' + + max_abs = 2.0 + tick_abs = 0.5 + max_anom = 2.0 + tick_anom = 1.0 + + deg2rad = pi/180. + lon_c = 50 + lat_c = -83 + radius = 10.5 + circle_bdry = -70+90 + + id = Dataset(paths[0]+file_tail, 'r') + data_tmp = id.variables['aice'][0,:350,:]*id.variables['hi'][0,:350,:] + lon_tmp = id.variables['TLON'][:350,:] + lat_tmp = id.variables['TLAT'][:350,:] + mask_tmp = id.variables['tmask'][:350,:] + id.close() + + lon = ma.empty([size(lon_tmp,0), size(lon_tmp,1)+1]) + lat = ma.empty([size(lat_tmp,0), size(lat_tmp,1)+1]) + mask = ma.empty([size(mask_tmp,0), size(mask_tmp,1)+1]) + data0 = ma.empty([size(data_tmp,0), size(data_tmp,1)+1]) + lon[:,:-1] = lon_tmp + lon[:,-1] = lon_tmp[:,0] + lat[:,:-1] = lat_tmp + lat[:,-1] = lat_tmp[:,0] + mask[:,:-1] = mask_tmp + mask[:,-1] = mask_tmp[:,0] + data0[:,:-1] = data_tmp + data0[:,-1] = data_tmp[:,0] + + land = ma.masked_where(mask==1, mask) + + x = -(lat+90)*cos(lon*deg2rad+pi/2) + y = (lat+90)*sin(lon*deg2rad+pi/2) + x_c = -(lat_c+90)*cos(lon_c*deg2rad+pi/2) + y_c = (lat_c+90)*sin(lon_c*deg2rad+pi/2) + x_reg, y_reg = meshgrid(linspace(-circle_bdry, circle_bdry, num=100), linspace(-circle_bdry, circle_bdry, num=100)) + land_circle = zeros(shape(x_reg)) + land_circle = ma.masked_where(sqrt((x_reg-x_c)**2 + (y_reg-y_c)**2) > radius, land_circle) + + fig = figure(figsize=(9,15)) + ax = fig.add_subplot(3, 2, 1, aspect='equal') + contourf(x, y, land, 1, colors=(('0.6', '0.6', '0.6'))) + contourf(x_reg, y_reg, land_circle, 1, colors=(('0.6', '0.6', '0.6'))) + img0 = pcolor(x, y, data0, vmin=0, vmax=max_abs, cmap='jet') #cmaps.viridis) + axis('off') + title(labels[0], fontsize=20) + cbaxes0 = fig.add_axes([0.025, 0.7, 0.02, 0.2]) + cbar0 = colorbar(img0, ticks=arange(0, max_abs+tick_abs, tick_abs), cax=cbaxes0, extend='max') + cbar0.ax.tick_params(labelsize=16) + + for sim in range(1, len(paths)): + id = Dataset(paths[sim]+file_tail, 'r') + data_tmp = id.variables['aice'][0,:350,:]*id.variables['hi'][0,:350,:] + id.close() + data = ma.empty([size(data_tmp,0), size(data_tmp,1)+1]) + data[:,:-1] = data_tmp + data[:,-1] = data_tmp[:,0] + data = data - data0 + ax = fig.add_subplot(3, 2, sim+1, aspect='equal') + contourf(x, y, land, 1, colors=(('0.6', '0.6', '0.6'))) + contourf(x_reg, y_reg, land_circle, 1, colors=(('0.6', '0.6', '0.6'))) + img = pcolor(x, y, data, vmin=-max_anom, vmax=max_anom, cmap='RdBu_r') + axis('off') + title(labels[sim], fontsize=20) + if sim == 3: + cbaxes = fig.add_axes([0.025, 0.4, 0.02, 0.2]) + cbar = colorbar(img, ticks=arange(-max_anom, max_anom+tick_anom, tick_anom), cax=cbaxes, extend='both') + cbar.ax.tick_params(labelsize=16) + + suptitle('Effective sea ice thickness (m) on 23 August', fontsize=28) + subplots_adjust(wspace=0.025,hspace=0.15) + + #fig.show() + fig.savefig('adv_thickness.png') + + +if __name__ == '__main__': + + adv_thickness() + diff --git a/adv_timeseries_volume.py b/adv_timeseries_volume.py index 253707f..95395f2 100644 --- a/adv_timeseries_volume.py +++ b/adv_timeseries_volume.py @@ -10,11 +10,11 @@ def adv_timeseries_volume (): # Number of output steps in the simulation (daily averages, 1 leap year) num_days = 366 # Paths to simulation directories (in order of decreasing sea ice volume) - paths = ['/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/c4_lowdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/a4_lowdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/u3_lowdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/c4_highdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/a4_highdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/u3_highdif/'] + paths = ['/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/c4_lowdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/a4_lowdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/u3_lowdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/c4_highdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/a4_highdif/', '/short/m68/kaa561/ROMS-CICE-MCT/tmproms/run/advection/u3limiters_lowdif/'] # Name of timeseries_seaice logfile in each directory logfile = 'seaice.log' # Abbreviations for each simulation to use in the legend - labels = ['C4_L', 'A4_L', 'U3_L', 'C4_H', 'A4_H', 'U3_H'] + labels = ['C4_L', 'A4_L', 'U3_L','C4_H', 'A4_H', 'U3_LIM'] # Colours for plotting each simulation colours = ['k', 'm', 'b', 'r', 'g', 'c'] diff --git a/aice_hi_seasonal.py b/aice_hi_seasonal.py index 30bb525..9e93373 100644 --- a/aice_hi_seasonal.py +++ b/aice_hi_seasonal.py @@ -2,242 +2,180 @@ from netCDF4 import Dataset, num2date from matplotlib.pyplot import * -# Creates a 4x2 plot of seasonally averaged sea ice concentration (top row) and -# thickness (bottom row) over the last year of simulation. -# Input: -# cice_file = path to CICE output file with 5-day averages, containing at least -# one complete Dec-Nov period (if there are multiple such periods, -# this script uses the last one) -# save = optional boolean to save the figure to a file, rather than displaying -# it on the screen -# fig_name = if save=True, path to the desired filename for the figure def aice_hi_seasonal (cice_file, save=False, fig_name=None): - # Starting and ending months (1-based) for each season start_month = [12, 3, 6, 9] end_month = [2, 5, 8, 11] - # Starting and ending days of the month (1-based) for each season - # Assume no leap years, we'll fix this later if we need start_day = [1, 1, 1, 1] end_day = [28, 31, 31, 30] - # Number of days in each season (again, ignore leap years for now) ndays_season = [90, 92, 92, 91] - # Season names for titles season_names = ['DJF', 'MAM', 'JJA', 'SON'] - # Degrees to radians conversion deg2rad = pi/180.0 - # Read CICE grid and time values id = Dataset(cice_file, 'r') - cice_lon_tmp = id.variables['TLON'][:-15,:] - cice_lat_tmp = id.variables['TLAT'][:-15,:] - # Wrap the periodic boundary by 1 cell - cice_lon = ma.empty([size(cice_lon_tmp,0), size(cice_lon_tmp,1)+1]) - cice_lat = ma.empty([size(cice_lat_tmp,0), size(cice_lat_tmp,1)+1]) - cice_lon[:,:-1] = cice_lon_tmp - cice_lon[:,-1] = cice_lon_tmp[:,0] - cice_lat[:,:-1] = cice_lat_tmp - cice_lat[:,-1] = cice_lat_tmp[:,0] + lon_tmp = id.variables['TLON'][:-15,:] + lat_tmp = id.variables['TLAT'][:-15,:] + lon = ma.empty([size(lon_tmp,0), size(lon_tmp,1)+1]) + lat = ma.empty([size(lat_tmp,0), size(lat_tmp,1)+1]) + lon[:,:-1] = lon_tmp + lon[:,-1] = lon_tmp[:,0] + lat[:,:-1] = lat_tmp + lat[:,-1] = lat_tmp[:,0] time_id = id.variables['time'] - # Get the year, month, and day (all 1-based) for each output step - # These are 5-day averages marked with the last day's date. - cice_time = num2date(time_id[:], units=time_id.units, calendar=time_id.calendar.lower()) - - # Loop backwards through time indices to find the last one we care about - # (which contains 30 November in its averaging period) - end_t = -1 # Missing value flag - for t in range(size(cice_time)-1, -1, -1): - if cice_time[t].month == end_month[-1] and cice_time[t].day == end_day[-1]: - end_t = t - break - if cice_time[t].month == start_month[0] and cice_time[t].day in range(start_day[0], start_day[0]+4): + time = num2date(time_id[:], units=time_id.units, calendar=time_id.calendar.lower()) + + end_t = -1 + for t in range(size(time)-1, -1, -1): + if time[t].month == start_month[0] and time[t].day in range(start_day[0], start_day[0]+5): end_t = t break - # Make sure we actually found it if end_t == -1: print 'Error: ' + cice_file + ' does not contain a complete Dec-Nov period' return - # Continue looping backwards to find the first time index we care about - # (which contains 1 December the previous year in its averaging period) - start_t = -1 # Missing value flag + start_t = -1 for t in range(end_t-60, -1, -1): - if cice_time[t].month == start_month[0] and cice_time[t].day in range(start_day[0], start_day[0]+5): + if time[t].month == start_month[0] and time[t].day in range(start_day[0]+1, start_day[0]+6): start_t = t break - # Make sure we actually found it if start_t == -1: print 'Error: ' + cice_file + ' does not contain a complete Dec-Nov period' return - # Check if end_t occurs on a leap year leap_year = False - if mod(cice_time[end_t].year, 4) == 0: - # Years divisible by 4 are leap years + if mod(time[end_t].year, 4) == 0: leap_year = True - if mod(cice_time[end_t].year, 100) == 0: - # Unless they're also divisible by 100, in which case they aren't - # leap years + if mod(time[end_t].year, 100) == 0: leap_year = False - if mod(cice_time[end_t].year, 400) == 0: - # Unless they're also divisible by 400, in which case they are - # leap years after all + if mod(time[end_t].year, 400) == 0: leap_year = True if leap_year: - # Update last day in February end_day[0] += 1 ndays_season[0] += 1 - # Initialise seasonal averages of CICE output - aice_tmp = ma.empty([4, size(cice_lon_tmp,0), size(cice_lon_tmp,1)]) + aice_tmp = ma.empty([4, size(lon_tmp,0), size(lon_tmp,1)]) aice_tmp[:,:,:] = 0.0 - hi_tmp = ma.empty([4, size(cice_lon_tmp,0), size(cice_lon_tmp,1)]) + hi_tmp = ma.empty([4, size(lon_tmp,0), size(lon_tmp,1)]) hi_tmp[:,:,:] = 0.0 - # Process one season at a time for season in range(4): - season_days = 0 # Number of days in season; this will be incremented + season_days = 0 next_season = mod(season+1, 4) - # Find starting timestep start_t_season = -1 for t in range(start_t, end_t+1): - if cice_time[t].month == start_month[season] and cice_time[t].day in range(start_day[season], start_day[season]+5): + if time[t].month == start_month[season] and time[t].day in range(start_day[season]+1, start_day[season]+6): start_t_season = t break - # Make sure we actually found it if start_t_season == -1: print 'Error: could not find starting timestep for season ' + season_names[season] return - # Find ending timestep end_t_season = -1 for t in range(start_t_season+1, end_t+1): - if cice_time[t].month == end_month[season] and cice_time[t].day == end_day[season]: + if time[t].month == start_month[next_season] and time[t].day in range(start_day[next_season], start_day[next_season]+5): end_t_season = t break - if cice_time[t].month == start_month[next_season] and cice_time[t].day in range(start_day[next_season], start_day[next_season]+4): - end_t_season = t - break - # Make sure we actually found it if end_t_season == -1: print 'Error: could not find ending timestep for season ' + season_names[season] return - # Figure out how many of the 5 days averaged in start_t_season are - # actually within this season - if cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 4: - # Starting day is in position 1 of 5; we care about all of them + if time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season] + 5: start_days = 5 - elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 3: - # Starting day is in position 2 of 5; we care about the last 4 + elif time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season] + 4: start_days = 4 - elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season]+ 2: - # Starting day is in position 3 of 5; we care about the last 3 + elif time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season]+ 3: start_days = 3 - elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 1: - # Starting day is in position 4 of 5; we care about the last 2 + elif time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season] + 2: start_days = 2 - elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season]: - # Starting day is in position 5 of 5; we care about the last 1 + elif time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season] + 1: start_days = 1 else: - print 'Error for season ' + season_names[season] + ': starting index is month ' + str(cice_time[start_t_season].month) + ', day ' + str(cice_time[start_t_season].day) + print 'Error for season ' + season_names[season] + ': starting index is month ' + str(time[start_t_season].month) + ', day ' + str(time[start_t_season].day) return - # Start accumulating data weighted by days aice_tmp[season,:,:] += id.variables['aice'][start_t_season,:-15,:]*start_days hi_tmp[season,:,:] += id.variables['hi'][start_t_season,:-15,:]*start_days season_days += start_days - # Beween start_t_season and end_t_season, we want all the days for t in range(start_t_season+1, end_t_season): aice_tmp[season,:,:] += id.variables['aice'][t,:-15,:]*5 hi_tmp[season,:,:] += id.variables['hi'][t,:-15,:]*5 season_days += 5 - # Figure out how many of the 5 days averaged in end_t_season are - # actually within this season - if cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 3: - # Ending day is in position 1 of 5; we care about the first 1 + if time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season] + 4: end_days = 1 - elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 2: - # Ending day is in position 2 of 5; we care about the first 2 + elif time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season] + 3: end_days = 2 - elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 1: - # Ending day is in position 3 of 5; we care about the first 3 + elif time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season] + 2: end_days = 3 - elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season]: - # Ending day is in position 4 of 5; we care about the first 4 + elif time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season] + 1: end_days = 4 - elif cice_time[end_t_season].month == end_month[season] and cice_time[end_t_season].day == end_day[season]: - # Ending day is in position 5 of 5; we care about all 5 + elif time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season]: end_days = 5 else: - print 'Error for season ' + season_names[season] + ': ending index is month ' + str(cice_time[end_t_season].month) + ', day ' + str(cice_time[end_t_season].day) + print 'Error for season ' + season_names[season] + ': ending index is month ' + str(time[end_t_season].month) + ', day ' + str(time[end_t_season].day) return aice_tmp[season,:,:] += id.variables['aice'][end_t_season,:-15,:]*end_days hi_tmp[season,:,:] += id.variables['hi'][end_t_season,:-15,:]*end_days season_days += end_days - # Check that we got the correct number of days if season_days != ndays_season[season]: print 'Error: found ' + str(season_days) + ' days instead of ' + str(ndays_season[season]) return - # Finished accumulating data, now convert from sum to average aice_tmp[season,:,:] /= season_days hi_tmp[season,:,:] /= season_days - # Finished reading all CICE data id.close() - # Wrap the periodic boundary aice = ma.empty([size(aice_tmp,0), size(aice_tmp,1), size(aice_tmp,2)+1]) - hi = ma.empty([size(hi_tmp,0), size(hi_tmp,1), size(hi_tmp,2)+1]) aice[:,:,:-1] = aice_tmp aice[:,:,-1] = aice_tmp[:,:,0] + hi = ma.empty([size(hi_tmp,0), size(hi_tmp,1), size(hi_tmp,2)+1]) hi[:,:,:-1] = hi_tmp - hi[:,:,-1] = hi_tmp[:,:,0] + hi[:,:,-1] = hi_tmp[:,:,0] + + x = -(lat+90)*cos(lon*deg2rad+pi/2) + y = (lat+90)*sin(lon*deg2rad+pi/2) - # Convert to spherical coordinates - cice_x = -(cice_lat+90)*cos(cice_lon*deg2rad+pi/2) - cice_y = (cice_lat+90)*sin(cice_lon*deg2rad+pi/2) + # Find boundaries for each side of plot based on extent of NSIDC grid + bdry1 = -35 + bdry2 = 35 + bdry3 = -33 + bdry4 = 37 # Set consistent colour levels lev1 = linspace(0, 1, num=50) - lev2 = linspace(0, 3, num=50) + lev2 = linspace(0, 2.5, num=50) # Plot fig = figure(figsize=(20,9)) # Loop over seasons for season in range(4): - # aice ax = fig.add_subplot(2, 4, season+1, aspect='equal') - img = contourf(cice_x, cice_y, aice[season,:,:], lev1, extend='both') + img = contourf(x, y, aice[season,:,:], lev1) + xlim([bdry1, bdry2]) + ylim([bdry3, bdry4]) + axis('off') if season == 0: - text(-39, 0, 'aice (%)', fontsize=21, ha='right') + text(-39, 0, 'aice (1)', fontsize=21, ha='right') title(season_names[season], fontsize=24) - xlim([-35, 35]) - ylim([-33, 37]) - axis('off') if season == 3: cbaxes1 = fig.add_axes([0.92, 0.55, 0.01, 0.3]) cbar1 = colorbar(img, ticks=arange(0,1+0.25,0.25), cax=cbaxes1) cbar1.ax.tick_params(labelsize=16) - # hi ax = fig.add_subplot(2, 4, season+5, aspect='equal') - img = contourf(cice_x, cice_y, hi[season,:,:], lev2, extend='both') + img = contourf(x, y, hi[season,:,:], lev2, extend='both') + xlim([bdry1, bdry2]) + ylim([bdry3, bdry4]) + axis('off') if season == 0: text(-39, 0, 'hi (m)', fontsize=21, ha='right') - xlim([-35, 35]) - ylim([-33, 37]) - axis('off') if season == 3: cbaxes2 = fig.add_axes([0.92, 0.15, 0.01, 0.3]) - cbar2 = colorbar(img, ticks=arange(0,3+1,1), cax=cbaxes2) + cbar2 = colorbar(img, ticks=arange(0,2.5+0.5,0.5), cax=cbaxes2) cbar2.ax.tick_params(labelsize=16) - # Decrease space between plots subplots_adjust(wspace=0.025,hspace=0.025) # Finished diff --git a/cice_vectorplot.py b/cice_vectorplot.py index e86e1e9..ae37e6e 100644 --- a/cice_vectorplot.py +++ b/cice_vectorplot.py @@ -56,30 +56,28 @@ def cice_vectorplot (file_path, tstep, xname, yname, cmax=None, save=False, fig_ speed = sqrt(u**2 + v**2) # Calculate X and Y coordinates for plotting circumpolar projection - X = -(lat+90)*cos(lon*deg2rad+pi/2) - Y = (lat+90)*sin(lon*deg2rad+pi/2) + x = -(lat+90)*cos(lon*deg2rad+pi/2) + y = (lat+90)*sin(lon*deg2rad+pi/2) - # Calculate vector components in spherical coordinate space - # (just differentiate and rearrange spherical coordinate transformation) - dlon_dt = u/(r*cos(lat*deg2rad)*deg2rad) - dlat_dt = v/(r*deg2rad) - dX_dt = -dlat_dt*cos(lon*deg2rad+pi/2) + (lat+90)*sin(lon*deg2rad+pi/2)*dlon_dt*deg2rad - dY_dt = dlat_dt*sin(lon*deg2rad+pi/2) + (lat+90)*cos(lon*deg2rad+pi/2)*dlon_dt*deg2rad + theta = arctan2(v, u) + theta_circ = theta - lon*deg2rad + u_circ = speed*cos(theta_circ) + v_circ = speed*sin(theta_circ) # Average X, Y, dX_dt, and dY_dt over block x block intervals # Calculate number of blocks - size0 = int(ceil(size(X,0)/float(block))) - size1 = int(ceil((size(X,1)-1)/float(block))) + size0 = int(ceil(size(x,0)/float(block))) + size1 = int(ceil((size(x,1)-1)/float(block))) # Set up arrays for averaged fields - X_block = ma.empty([size0, size1]) - Y_block = ma.empty([size0, size1]) - dX_dt_block = ma.empty([size0, size1]) - dY_dt_block = ma.empty([size0, size1]) + x_block = ma.empty([size0, size1]) + y_block = ma.empty([size0, size1]) + u_circ_block = ma.empty([size0, size1]) + v_circ_block = ma.empty([size0, size1]) # Set up arrays containing boundary indices - posn0 = range(0, size(X,0), block) - posn0.append(size(X,0)) - posn1 = range(0, size(X,1), block) - posn1.append(size(X,1)) + posn0 = range(0, size(x,0), block) + posn0.append(size(x,0)) + posn1 = range(0, size(x,1), block) + posn1.append(size(x,1)) # Double loop to average each block (can't find a more efficient way to do # this) for j in range(size0): @@ -88,10 +86,10 @@ def cice_vectorplot (file_path, tstep, xname, yname, cmax=None, save=False, fig_ end0 = posn0[j+1] start1 = posn1[i] end1 = posn1[i+1] - X_block[j,i] = mean(X[start0:end0, start1:end1]) - Y_block[j,i] = mean(Y[start0:end0, start1:end1]) - dX_dt_block[j,i] = mean(dX_dt[start0:end0, start1:end1]) - dY_dt_block[j,i] = mean(dY_dt[start0:end0, start1:end1]) + x_block[j,i] = mean(x[start0:end0, start1:end1]) + y_block[j,i] = mean(y[start0:end0, start1:end1]) + u_circ_block[j,i] = mean(u_circ[start0:end0, start1:end1]) + v_circ_block[j,i] = mean(v_circ[start0:end0, start1:end1]) # Set up colour scale levels if cmax is None: @@ -104,11 +102,11 @@ def cice_vectorplot (file_path, tstep, xname, yname, cmax=None, save=False, fig_ fig.add_subplot(1,1,1, aspect='equal') # Contour speed values at every point # Use pastel colour map so overlaid vectors will show up - contourf(X, Y, speed, lev, cmap='Paired', extend='both') + contourf(x, y, speed, lev, cmap='Paired', extend='both') cbar = colorbar() cbar.ax.tick_params(labelsize=20) # Add vectors for each block - quiver(X_block, Y_block, dX_dt_block, dY_dt_block, color='black') + quiver(x_block, y_block, u_circ_block, v_circ_block, color='black') title(xname + ', ' + yname, fontsize=30) axis('off') diff --git a/colormaps.py b/colormaps.py new file mode 100644 index 0000000..d1f74a4 --- /dev/null +++ b/colormaps.py @@ -0,0 +1,1058 @@ +# New matplotlib colormaps by Nathaniel J. Smith, Stefan van der Walt, +# and (in the case of viridis) Eric Firing. +# +# This file and the colormaps in it are released under the CC0 license / +# public domain dedication. We would appreciate credit if you use or +# redistribute these colormaps, but do not impose any legal restrictions. +# +# To the extent possible under law, the persons who associated CC0 with +# mpl-colormaps have waived all copyright and related or neighboring rights +# to mpl-colormaps. +# +# You should have received a copy of the CC0 legalcode along with this +# work. If not, see . + +__all__ = ['magma', 'inferno', 'plasma', 'viridis'] + +_magma_data = [[0.001462, 0.000466, 0.013866], + [0.002258, 0.001295, 0.018331], + [0.003279, 0.002305, 0.023708], + [0.004512, 0.003490, 0.029965], + [0.005950, 0.004843, 0.037130], + [0.007588, 0.006356, 0.044973], + [0.009426, 0.008022, 0.052844], + [0.011465, 0.009828, 0.060750], + [0.013708, 0.011771, 0.068667], + [0.016156, 0.013840, 0.076603], + [0.018815, 0.016026, 0.084584], + [0.021692, 0.018320, 0.092610], + [0.024792, 0.020715, 0.100676], + [0.028123, 0.023201, 0.108787], + [0.031696, 0.025765, 0.116965], + [0.035520, 0.028397, 0.125209], + [0.039608, 0.031090, 0.133515], + [0.043830, 0.033830, 0.141886], + [0.048062, 0.036607, 0.150327], + [0.052320, 0.039407, 0.158841], + [0.056615, 0.042160, 0.167446], + [0.060949, 0.044794, 0.176129], + [0.065330, 0.047318, 0.184892], + [0.069764, 0.049726, 0.193735], + [0.074257, 0.052017, 0.202660], + [0.078815, 0.054184, 0.211667], + [0.083446, 0.056225, 0.220755], + [0.088155, 0.058133, 0.229922], + [0.092949, 0.059904, 0.239164], + [0.097833, 0.061531, 0.248477], + [0.102815, 0.063010, 0.257854], + [0.107899, 0.064335, 0.267289], + [0.113094, 0.065492, 0.276784], + [0.118405, 0.066479, 0.286321], + [0.123833, 0.067295, 0.295879], + [0.129380, 0.067935, 0.305443], + [0.135053, 0.068391, 0.315000], + [0.140858, 0.068654, 0.324538], + [0.146785, 0.068738, 0.334011], + [0.152839, 0.068637, 0.343404], + [0.159018, 0.068354, 0.352688], + [0.165308, 0.067911, 0.361816], + [0.171713, 0.067305, 0.370771], + [0.178212, 0.066576, 0.379497], + [0.184801, 0.065732, 0.387973], + [0.191460, 0.064818, 0.396152], + [0.198177, 0.063862, 0.404009], + [0.204935, 0.062907, 0.411514], + [0.211718, 0.061992, 0.418647], + [0.218512, 0.061158, 0.425392], + [0.225302, 0.060445, 0.431742], + [0.232077, 0.059889, 0.437695], + [0.238826, 0.059517, 0.443256], + [0.245543, 0.059352, 0.448436], + [0.252220, 0.059415, 0.453248], + [0.258857, 0.059706, 0.457710], + [0.265447, 0.060237, 0.461840], + [0.271994, 0.060994, 0.465660], + [0.278493, 0.061978, 0.469190], + [0.284951, 0.063168, 0.472451], + [0.291366, 0.064553, 0.475462], + [0.297740, 0.066117, 0.478243], + [0.304081, 0.067835, 0.480812], + [0.310382, 0.069702, 0.483186], + [0.316654, 0.071690, 0.485380], + [0.322899, 0.073782, 0.487408], + [0.329114, 0.075972, 0.489287], + [0.335308, 0.078236, 0.491024], + [0.341482, 0.080564, 0.492631], + [0.347636, 0.082946, 0.494121], + [0.353773, 0.085373, 0.495501], + [0.359898, 0.087831, 0.496778], + [0.366012, 0.090314, 0.497960], + [0.372116, 0.092816, 0.499053], + [0.378211, 0.095332, 0.500067], + [0.384299, 0.097855, 0.501002], + [0.390384, 0.100379, 0.501864], + [0.396467, 0.102902, 0.502658], + [0.402548, 0.105420, 0.503386], + [0.408629, 0.107930, 0.504052], + [0.414709, 0.110431, 0.504662], + [0.420791, 0.112920, 0.505215], + [0.426877, 0.115395, 0.505714], + [0.432967, 0.117855, 0.506160], + [0.439062, 0.120298, 0.506555], + [0.445163, 0.122724, 0.506901], + [0.451271, 0.125132, 0.507198], + [0.457386, 0.127522, 0.507448], + [0.463508, 0.129893, 0.507652], + [0.469640, 0.132245, 0.507809], + [0.475780, 0.134577, 0.507921], + [0.481929, 0.136891, 0.507989], + [0.488088, 0.139186, 0.508011], + [0.494258, 0.141462, 0.507988], + [0.500438, 0.143719, 0.507920], + [0.506629, 0.145958, 0.507806], + [0.512831, 0.148179, 0.507648], + [0.519045, 0.150383, 0.507443], + [0.525270, 0.152569, 0.507192], + [0.531507, 0.154739, 0.506895], + [0.537755, 0.156894, 0.506551], + [0.544015, 0.159033, 0.506159], + [0.550287, 0.161158, 0.505719], + [0.556571, 0.163269, 0.505230], + [0.562866, 0.165368, 0.504692], + [0.569172, 0.167454, 0.504105], + [0.575490, 0.169530, 0.503466], + [0.581819, 0.171596, 0.502777], + [0.588158, 0.173652, 0.502035], + [0.594508, 0.175701, 0.501241], + [0.600868, 0.177743, 0.500394], + [0.607238, 0.179779, 0.499492], + [0.613617, 0.181811, 0.498536], + [0.620005, 0.183840, 0.497524], + [0.626401, 0.185867, 0.496456], + [0.632805, 0.187893, 0.495332], + [0.639216, 0.189921, 0.494150], + [0.645633, 0.191952, 0.492910], + [0.652056, 0.193986, 0.491611], + [0.658483, 0.196027, 0.490253], + [0.664915, 0.198075, 0.488836], + [0.671349, 0.200133, 0.487358], + [0.677786, 0.202203, 0.485819], + [0.684224, 0.204286, 0.484219], + [0.690661, 0.206384, 0.482558], + [0.697098, 0.208501, 0.480835], + [0.703532, 0.210638, 0.479049], + [0.709962, 0.212797, 0.477201], + [0.716387, 0.214982, 0.475290], + [0.722805, 0.217194, 0.473316], + [0.729216, 0.219437, 0.471279], + [0.735616, 0.221713, 0.469180], + [0.742004, 0.224025, 0.467018], + [0.748378, 0.226377, 0.464794], + [0.754737, 0.228772, 0.462509], + [0.761077, 0.231214, 0.460162], + [0.767398, 0.233705, 0.457755], + [0.773695, 0.236249, 0.455289], + [0.779968, 0.238851, 0.452765], + [0.786212, 0.241514, 0.450184], + [0.792427, 0.244242, 0.447543], + [0.798608, 0.247040, 0.444848], + [0.804752, 0.249911, 0.442102], + [0.810855, 0.252861, 0.439305], + [0.816914, 0.255895, 0.436461], + [0.822926, 0.259016, 0.433573], + [0.828886, 0.262229, 0.430644], + [0.834791, 0.265540, 0.427671], + [0.840636, 0.268953, 0.424666], + [0.846416, 0.272473, 0.421631], + [0.852126, 0.276106, 0.418573], + [0.857763, 0.279857, 0.415496], + [0.863320, 0.283729, 0.412403], + [0.868793, 0.287728, 0.409303], + [0.874176, 0.291859, 0.406205], + [0.879464, 0.296125, 0.403118], + [0.884651, 0.300530, 0.400047], + [0.889731, 0.305079, 0.397002], + [0.894700, 0.309773, 0.393995], + [0.899552, 0.314616, 0.391037], + [0.904281, 0.319610, 0.388137], + [0.908884, 0.324755, 0.385308], + [0.913354, 0.330052, 0.382563], + [0.917689, 0.335500, 0.379915], + [0.921884, 0.341098, 0.377376], + [0.925937, 0.346844, 0.374959], + [0.929845, 0.352734, 0.372677], + [0.933606, 0.358764, 0.370541], + [0.937221, 0.364929, 0.368567], + [0.940687, 0.371224, 0.366762], + [0.944006, 0.377643, 0.365136], + [0.947180, 0.384178, 0.363701], + [0.950210, 0.390820, 0.362468], + [0.953099, 0.397563, 0.361438], + [0.955849, 0.404400, 0.360619], + [0.958464, 0.411324, 0.360014], + [0.960949, 0.418323, 0.359630], + [0.963310, 0.425390, 0.359469], + [0.965549, 0.432519, 0.359529], + [0.967671, 0.439703, 0.359810], + [0.969680, 0.446936, 0.360311], + [0.971582, 0.454210, 0.361030], + [0.973381, 0.461520, 0.361965], + [0.975082, 0.468861, 0.363111], + [0.976690, 0.476226, 0.364466], + [0.978210, 0.483612, 0.366025], + [0.979645, 0.491014, 0.367783], + [0.981000, 0.498428, 0.369734], + [0.982279, 0.505851, 0.371874], + [0.983485, 0.513280, 0.374198], + [0.984622, 0.520713, 0.376698], + [0.985693, 0.528148, 0.379371], + [0.986700, 0.535582, 0.382210], + [0.987646, 0.543015, 0.385210], + [0.988533, 0.550446, 0.388365], + [0.989363, 0.557873, 0.391671], + [0.990138, 0.565296, 0.395122], + [0.990871, 0.572706, 0.398714], + [0.991558, 0.580107, 0.402441], + [0.992196, 0.587502, 0.406299], + [0.992785, 0.594891, 0.410283], + [0.993326, 0.602275, 0.414390], + [0.993834, 0.609644, 0.418613], + [0.994309, 0.616999, 0.422950], + [0.994738, 0.624350, 0.427397], + [0.995122, 0.631696, 0.431951], + [0.995480, 0.639027, 0.436607], + [0.995810, 0.646344, 0.441361], + [0.996096, 0.653659, 0.446213], + [0.996341, 0.660969, 0.451160], + [0.996580, 0.668256, 0.456192], + [0.996775, 0.675541, 0.461314], + [0.996925, 0.682828, 0.466526], + [0.997077, 0.690088, 0.471811], + [0.997186, 0.697349, 0.477182], + [0.997254, 0.704611, 0.482635], + [0.997325, 0.711848, 0.488154], + [0.997351, 0.719089, 0.493755], + [0.997351, 0.726324, 0.499428], + [0.997341, 0.733545, 0.505167], + [0.997285, 0.740772, 0.510983], + [0.997228, 0.747981, 0.516859], + [0.997138, 0.755190, 0.522806], + [0.997019, 0.762398, 0.528821], + [0.996898, 0.769591, 0.534892], + [0.996727, 0.776795, 0.541039], + [0.996571, 0.783977, 0.547233], + [0.996369, 0.791167, 0.553499], + [0.996162, 0.798348, 0.559820], + [0.995932, 0.805527, 0.566202], + [0.995680, 0.812706, 0.572645], + [0.995424, 0.819875, 0.579140], + [0.995131, 0.827052, 0.585701], + [0.994851, 0.834213, 0.592307], + [0.994524, 0.841387, 0.598983], + [0.994222, 0.848540, 0.605696], + [0.993866, 0.855711, 0.612482], + [0.993545, 0.862859, 0.619299], + [0.993170, 0.870024, 0.626189], + [0.992831, 0.877168, 0.633109], + [0.992440, 0.884330, 0.640099], + [0.992089, 0.891470, 0.647116], + [0.991688, 0.898627, 0.654202], + [0.991332, 0.905763, 0.661309], + [0.990930, 0.912915, 0.668481], + [0.990570, 0.920049, 0.675675], + [0.990175, 0.927196, 0.682926], + [0.989815, 0.934329, 0.690198], + [0.989434, 0.941470, 0.697519], + [0.989077, 0.948604, 0.704863], + [0.988717, 0.955742, 0.712242], + [0.988367, 0.962878, 0.719649], + [0.988033, 0.970012, 0.727077], + [0.987691, 0.977154, 0.734536], + [0.987387, 0.984288, 0.742002], + [0.987053, 0.991438, 0.749504]] + +_inferno_data = [[0.001462, 0.000466, 0.013866], + [0.002267, 0.001270, 0.018570], + [0.003299, 0.002249, 0.024239], + [0.004547, 0.003392, 0.030909], + [0.006006, 0.004692, 0.038558], + [0.007676, 0.006136, 0.046836], + [0.009561, 0.007713, 0.055143], + [0.011663, 0.009417, 0.063460], + [0.013995, 0.011225, 0.071862], + [0.016561, 0.013136, 0.080282], + [0.019373, 0.015133, 0.088767], + [0.022447, 0.017199, 0.097327], + [0.025793, 0.019331, 0.105930], + [0.029432, 0.021503, 0.114621], + [0.033385, 0.023702, 0.123397], + [0.037668, 0.025921, 0.132232], + [0.042253, 0.028139, 0.141141], + [0.046915, 0.030324, 0.150164], + [0.051644, 0.032474, 0.159254], + [0.056449, 0.034569, 0.168414], + [0.061340, 0.036590, 0.177642], + [0.066331, 0.038504, 0.186962], + [0.071429, 0.040294, 0.196354], + [0.076637, 0.041905, 0.205799], + [0.081962, 0.043328, 0.215289], + [0.087411, 0.044556, 0.224813], + [0.092990, 0.045583, 0.234358], + [0.098702, 0.046402, 0.243904], + [0.104551, 0.047008, 0.253430], + [0.110536, 0.047399, 0.262912], + [0.116656, 0.047574, 0.272321], + [0.122908, 0.047536, 0.281624], + [0.129285, 0.047293, 0.290788], + [0.135778, 0.046856, 0.299776], + [0.142378, 0.046242, 0.308553], + [0.149073, 0.045468, 0.317085], + [0.155850, 0.044559, 0.325338], + [0.162689, 0.043554, 0.333277], + [0.169575, 0.042489, 0.340874], + [0.176493, 0.041402, 0.348111], + [0.183429, 0.040329, 0.354971], + [0.190367, 0.039309, 0.361447], + [0.197297, 0.038400, 0.367535], + [0.204209, 0.037632, 0.373238], + [0.211095, 0.037030, 0.378563], + [0.217949, 0.036615, 0.383522], + [0.224763, 0.036405, 0.388129], + [0.231538, 0.036405, 0.392400], + [0.238273, 0.036621, 0.396353], + [0.244967, 0.037055, 0.400007], + [0.251620, 0.037705, 0.403378], + [0.258234, 0.038571, 0.406485], + [0.264810, 0.039647, 0.409345], + [0.271347, 0.040922, 0.411976], + [0.277850, 0.042353, 0.414392], + [0.284321, 0.043933, 0.416608], + [0.290763, 0.045644, 0.418637], + [0.297178, 0.047470, 0.420491], + [0.303568, 0.049396, 0.422182], + [0.309935, 0.051407, 0.423721], + [0.316282, 0.053490, 0.425116], + [0.322610, 0.055634, 0.426377], + [0.328921, 0.057827, 0.427511], + [0.335217, 0.060060, 0.428524], + [0.341500, 0.062325, 0.429425], + [0.347771, 0.064616, 0.430217], + [0.354032, 0.066925, 0.430906], + [0.360284, 0.069247, 0.431497], + [0.366529, 0.071579, 0.431994], + [0.372768, 0.073915, 0.432400], + [0.379001, 0.076253, 0.432719], + [0.385228, 0.078591, 0.432955], + [0.391453, 0.080927, 0.433109], + [0.397674, 0.083257, 0.433183], + [0.403894, 0.085580, 0.433179], + [0.410113, 0.087896, 0.433098], + [0.416331, 0.090203, 0.432943], + [0.422549, 0.092501, 0.432714], + [0.428768, 0.094790, 0.432412], + [0.434987, 0.097069, 0.432039], + [0.441207, 0.099338, 0.431594], + [0.447428, 0.101597, 0.431080], + [0.453651, 0.103848, 0.430498], + [0.459875, 0.106089, 0.429846], + [0.466100, 0.108322, 0.429125], + [0.472328, 0.110547, 0.428334], + [0.478558, 0.112764, 0.427475], + [0.484789, 0.114974, 0.426548], + [0.491022, 0.117179, 0.425552], + [0.497257, 0.119379, 0.424488], + [0.503493, 0.121575, 0.423356], + [0.509730, 0.123769, 0.422156], + [0.515967, 0.125960, 0.420887], + [0.522206, 0.128150, 0.419549], + [0.528444, 0.130341, 0.418142], + [0.534683, 0.132534, 0.416667], + [0.540920, 0.134729, 0.415123], + [0.547157, 0.136929, 0.413511], + [0.553392, 0.139134, 0.411829], + [0.559624, 0.141346, 0.410078], + [0.565854, 0.143567, 0.408258], + [0.572081, 0.145797, 0.406369], + [0.578304, 0.148039, 0.404411], + [0.584521, 0.150294, 0.402385], + [0.590734, 0.152563, 0.400290], + [0.596940, 0.154848, 0.398125], + [0.603139, 0.157151, 0.395891], + [0.609330, 0.159474, 0.393589], + [0.615513, 0.161817, 0.391219], + [0.621685, 0.164184, 0.388781], + [0.627847, 0.166575, 0.386276], + [0.633998, 0.168992, 0.383704], + [0.640135, 0.171438, 0.381065], + [0.646260, 0.173914, 0.378359], + [0.652369, 0.176421, 0.375586], + [0.658463, 0.178962, 0.372748], + [0.664540, 0.181539, 0.369846], + [0.670599, 0.184153, 0.366879], + [0.676638, 0.186807, 0.363849], + [0.682656, 0.189501, 0.360757], + [0.688653, 0.192239, 0.357603], + [0.694627, 0.195021, 0.354388], + [0.700576, 0.197851, 0.351113], + [0.706500, 0.200728, 0.347777], + [0.712396, 0.203656, 0.344383], + [0.718264, 0.206636, 0.340931], + [0.724103, 0.209670, 0.337424], + [0.729909, 0.212759, 0.333861], + [0.735683, 0.215906, 0.330245], + [0.741423, 0.219112, 0.326576], + [0.747127, 0.222378, 0.322856], + [0.752794, 0.225706, 0.319085], + [0.758422, 0.229097, 0.315266], + [0.764010, 0.232554, 0.311399], + [0.769556, 0.236077, 0.307485], + [0.775059, 0.239667, 0.303526], + [0.780517, 0.243327, 0.299523], + [0.785929, 0.247056, 0.295477], + [0.791293, 0.250856, 0.291390], + [0.796607, 0.254728, 0.287264], + [0.801871, 0.258674, 0.283099], + [0.807082, 0.262692, 0.278898], + [0.812239, 0.266786, 0.274661], + [0.817341, 0.270954, 0.270390], + [0.822386, 0.275197, 0.266085], + [0.827372, 0.279517, 0.261750], + [0.832299, 0.283913, 0.257383], + [0.837165, 0.288385, 0.252988], + [0.841969, 0.292933, 0.248564], + [0.846709, 0.297559, 0.244113], + [0.851384, 0.302260, 0.239636], + [0.855992, 0.307038, 0.235133], + [0.860533, 0.311892, 0.230606], + [0.865006, 0.316822, 0.226055], + [0.869409, 0.321827, 0.221482], + [0.873741, 0.326906, 0.216886], + [0.878001, 0.332060, 0.212268], + [0.882188, 0.337287, 0.207628], + [0.886302, 0.342586, 0.202968], + [0.890341, 0.347957, 0.198286], + [0.894305, 0.353399, 0.193584], + [0.898192, 0.358911, 0.188860], + [0.902003, 0.364492, 0.184116], + [0.905735, 0.370140, 0.179350], + [0.909390, 0.375856, 0.174563], + [0.912966, 0.381636, 0.169755], + [0.916462, 0.387481, 0.164924], + [0.919879, 0.393389, 0.160070], + [0.923215, 0.399359, 0.155193], + [0.926470, 0.405389, 0.150292], + [0.929644, 0.411479, 0.145367], + [0.932737, 0.417627, 0.140417], + [0.935747, 0.423831, 0.135440], + [0.938675, 0.430091, 0.130438], + [0.941521, 0.436405, 0.125409], + [0.944285, 0.442772, 0.120354], + [0.946965, 0.449191, 0.115272], + [0.949562, 0.455660, 0.110164], + [0.952075, 0.462178, 0.105031], + [0.954506, 0.468744, 0.099874], + [0.956852, 0.475356, 0.094695], + [0.959114, 0.482014, 0.089499], + [0.961293, 0.488716, 0.084289], + [0.963387, 0.495462, 0.079073], + [0.965397, 0.502249, 0.073859], + [0.967322, 0.509078, 0.068659], + [0.969163, 0.515946, 0.063488], + [0.970919, 0.522853, 0.058367], + [0.972590, 0.529798, 0.053324], + [0.974176, 0.536780, 0.048392], + [0.975677, 0.543798, 0.043618], + [0.977092, 0.550850, 0.039050], + [0.978422, 0.557937, 0.034931], + [0.979666, 0.565057, 0.031409], + [0.980824, 0.572209, 0.028508], + [0.981895, 0.579392, 0.026250], + [0.982881, 0.586606, 0.024661], + [0.983779, 0.593849, 0.023770], + [0.984591, 0.601122, 0.023606], + [0.985315, 0.608422, 0.024202], + [0.985952, 0.615750, 0.025592], + [0.986502, 0.623105, 0.027814], + [0.986964, 0.630485, 0.030908], + [0.987337, 0.637890, 0.034916], + [0.987622, 0.645320, 0.039886], + [0.987819, 0.652773, 0.045581], + [0.987926, 0.660250, 0.051750], + [0.987945, 0.667748, 0.058329], + [0.987874, 0.675267, 0.065257], + [0.987714, 0.682807, 0.072489], + [0.987464, 0.690366, 0.079990], + [0.987124, 0.697944, 0.087731], + [0.986694, 0.705540, 0.095694], + [0.986175, 0.713153, 0.103863], + [0.985566, 0.720782, 0.112229], + [0.984865, 0.728427, 0.120785], + [0.984075, 0.736087, 0.129527], + [0.983196, 0.743758, 0.138453], + [0.982228, 0.751442, 0.147565], + [0.981173, 0.759135, 0.156863], + [0.980032, 0.766837, 0.166353], + [0.978806, 0.774545, 0.176037], + [0.977497, 0.782258, 0.185923], + [0.976108, 0.789974, 0.196018], + [0.974638, 0.797692, 0.206332], + [0.973088, 0.805409, 0.216877], + [0.971468, 0.813122, 0.227658], + [0.969783, 0.820825, 0.238686], + [0.968041, 0.828515, 0.249972], + [0.966243, 0.836191, 0.261534], + [0.964394, 0.843848, 0.273391], + [0.962517, 0.851476, 0.285546], + [0.960626, 0.859069, 0.298010], + [0.958720, 0.866624, 0.310820], + [0.956834, 0.874129, 0.323974], + [0.954997, 0.881569, 0.337475], + [0.953215, 0.888942, 0.351369], + [0.951546, 0.896226, 0.365627], + [0.950018, 0.903409, 0.380271], + [0.948683, 0.910473, 0.395289], + [0.947594, 0.917399, 0.410665], + [0.946809, 0.924168, 0.426373], + [0.946392, 0.930761, 0.442367], + [0.946403, 0.937159, 0.458592], + [0.946903, 0.943348, 0.474970], + [0.947937, 0.949318, 0.491426], + [0.949545, 0.955063, 0.507860], + [0.951740, 0.960587, 0.524203], + [0.954529, 0.965896, 0.540361], + [0.957896, 0.971003, 0.556275], + [0.961812, 0.975924, 0.571925], + [0.966249, 0.980678, 0.587206], + [0.971162, 0.985282, 0.602154], + [0.976511, 0.989753, 0.616760], + [0.982257, 0.994109, 0.631017], + [0.988362, 0.998364, 0.644924]] + +_plasma_data = [[0.050383, 0.029803, 0.527975], + [0.063536, 0.028426, 0.533124], + [0.075353, 0.027206, 0.538007], + [0.086222, 0.026125, 0.542658], + [0.096379, 0.025165, 0.547103], + [0.105980, 0.024309, 0.551368], + [0.115124, 0.023556, 0.555468], + [0.123903, 0.022878, 0.559423], + [0.132381, 0.022258, 0.563250], + [0.140603, 0.021687, 0.566959], + [0.148607, 0.021154, 0.570562], + [0.156421, 0.020651, 0.574065], + [0.164070, 0.020171, 0.577478], + [0.171574, 0.019706, 0.580806], + [0.178950, 0.019252, 0.584054], + [0.186213, 0.018803, 0.587228], + [0.193374, 0.018354, 0.590330], + [0.200445, 0.017902, 0.593364], + [0.207435, 0.017442, 0.596333], + [0.214350, 0.016973, 0.599239], + [0.221197, 0.016497, 0.602083], + [0.227983, 0.016007, 0.604867], + [0.234715, 0.015502, 0.607592], + [0.241396, 0.014979, 0.610259], + [0.248032, 0.014439, 0.612868], + [0.254627, 0.013882, 0.615419], + [0.261183, 0.013308, 0.617911], + [0.267703, 0.012716, 0.620346], + [0.274191, 0.012109, 0.622722], + [0.280648, 0.011488, 0.625038], + [0.287076, 0.010855, 0.627295], + [0.293478, 0.010213, 0.629490], + [0.299855, 0.009561, 0.631624], + [0.306210, 0.008902, 0.633694], + [0.312543, 0.008239, 0.635700], + [0.318856, 0.007576, 0.637640], + [0.325150, 0.006915, 0.639512], + [0.331426, 0.006261, 0.641316], + [0.337683, 0.005618, 0.643049], + [0.343925, 0.004991, 0.644710], + [0.350150, 0.004382, 0.646298], + [0.356359, 0.003798, 0.647810], + [0.362553, 0.003243, 0.649245], + [0.368733, 0.002724, 0.650601], + [0.374897, 0.002245, 0.651876], + [0.381047, 0.001814, 0.653068], + [0.387183, 0.001434, 0.654177], + [0.393304, 0.001114, 0.655199], + [0.399411, 0.000859, 0.656133], + [0.405503, 0.000678, 0.656977], + [0.411580, 0.000577, 0.657730], + [0.417642, 0.000564, 0.658390], + [0.423689, 0.000646, 0.658956], + [0.429719, 0.000831, 0.659425], + [0.435734, 0.001127, 0.659797], + [0.441732, 0.001540, 0.660069], + [0.447714, 0.002080, 0.660240], + [0.453677, 0.002755, 0.660310], + [0.459623, 0.003574, 0.660277], + [0.465550, 0.004545, 0.660139], + [0.471457, 0.005678, 0.659897], + [0.477344, 0.006980, 0.659549], + [0.483210, 0.008460, 0.659095], + [0.489055, 0.010127, 0.658534], + [0.494877, 0.011990, 0.657865], + [0.500678, 0.014055, 0.657088], + [0.506454, 0.016333, 0.656202], + [0.512206, 0.018833, 0.655209], + [0.517933, 0.021563, 0.654109], + [0.523633, 0.024532, 0.652901], + [0.529306, 0.027747, 0.651586], + [0.534952, 0.031217, 0.650165], + [0.540570, 0.034950, 0.648640], + [0.546157, 0.038954, 0.647010], + [0.551715, 0.043136, 0.645277], + [0.557243, 0.047331, 0.643443], + [0.562738, 0.051545, 0.641509], + [0.568201, 0.055778, 0.639477], + [0.573632, 0.060028, 0.637349], + [0.579029, 0.064296, 0.635126], + [0.584391, 0.068579, 0.632812], + [0.589719, 0.072878, 0.630408], + [0.595011, 0.077190, 0.627917], + [0.600266, 0.081516, 0.625342], + [0.605485, 0.085854, 0.622686], + [0.610667, 0.090204, 0.619951], + [0.615812, 0.094564, 0.617140], + [0.620919, 0.098934, 0.614257], + [0.625987, 0.103312, 0.611305], + [0.631017, 0.107699, 0.608287], + [0.636008, 0.112092, 0.605205], + [0.640959, 0.116492, 0.602065], + [0.645872, 0.120898, 0.598867], + [0.650746, 0.125309, 0.595617], + [0.655580, 0.129725, 0.592317], + [0.660374, 0.134144, 0.588971], + [0.665129, 0.138566, 0.585582], + [0.669845, 0.142992, 0.582154], + [0.674522, 0.147419, 0.578688], + [0.679160, 0.151848, 0.575189], + [0.683758, 0.156278, 0.571660], + [0.688318, 0.160709, 0.568103], + [0.692840, 0.165141, 0.564522], + [0.697324, 0.169573, 0.560919], + [0.701769, 0.174005, 0.557296], + [0.706178, 0.178437, 0.553657], + [0.710549, 0.182868, 0.550004], + [0.714883, 0.187299, 0.546338], + [0.719181, 0.191729, 0.542663], + [0.723444, 0.196158, 0.538981], + [0.727670, 0.200586, 0.535293], + [0.731862, 0.205013, 0.531601], + [0.736019, 0.209439, 0.527908], + [0.740143, 0.213864, 0.524216], + [0.744232, 0.218288, 0.520524], + [0.748289, 0.222711, 0.516834], + [0.752312, 0.227133, 0.513149], + [0.756304, 0.231555, 0.509468], + [0.760264, 0.235976, 0.505794], + [0.764193, 0.240396, 0.502126], + [0.768090, 0.244817, 0.498465], + [0.771958, 0.249237, 0.494813], + [0.775796, 0.253658, 0.491171], + [0.779604, 0.258078, 0.487539], + [0.783383, 0.262500, 0.483918], + [0.787133, 0.266922, 0.480307], + [0.790855, 0.271345, 0.476706], + [0.794549, 0.275770, 0.473117], + [0.798216, 0.280197, 0.469538], + [0.801855, 0.284626, 0.465971], + [0.805467, 0.289057, 0.462415], + [0.809052, 0.293491, 0.458870], + [0.812612, 0.297928, 0.455338], + [0.816144, 0.302368, 0.451816], + [0.819651, 0.306812, 0.448306], + [0.823132, 0.311261, 0.444806], + [0.826588, 0.315714, 0.441316], + [0.830018, 0.320172, 0.437836], + [0.833422, 0.324635, 0.434366], + [0.836801, 0.329105, 0.430905], + [0.840155, 0.333580, 0.427455], + [0.843484, 0.338062, 0.424013], + [0.846788, 0.342551, 0.420579], + [0.850066, 0.347048, 0.417153], + [0.853319, 0.351553, 0.413734], + [0.856547, 0.356066, 0.410322], + [0.859750, 0.360588, 0.406917], + [0.862927, 0.365119, 0.403519], + [0.866078, 0.369660, 0.400126], + [0.869203, 0.374212, 0.396738], + [0.872303, 0.378774, 0.393355], + [0.875376, 0.383347, 0.389976], + [0.878423, 0.387932, 0.386600], + [0.881443, 0.392529, 0.383229], + [0.884436, 0.397139, 0.379860], + [0.887402, 0.401762, 0.376494], + [0.890340, 0.406398, 0.373130], + [0.893250, 0.411048, 0.369768], + [0.896131, 0.415712, 0.366407], + [0.898984, 0.420392, 0.363047], + [0.901807, 0.425087, 0.359688], + [0.904601, 0.429797, 0.356329], + [0.907365, 0.434524, 0.352970], + [0.910098, 0.439268, 0.349610], + [0.912800, 0.444029, 0.346251], + [0.915471, 0.448807, 0.342890], + [0.918109, 0.453603, 0.339529], + [0.920714, 0.458417, 0.336166], + [0.923287, 0.463251, 0.332801], + [0.925825, 0.468103, 0.329435], + [0.928329, 0.472975, 0.326067], + [0.930798, 0.477867, 0.322697], + [0.933232, 0.482780, 0.319325], + [0.935630, 0.487712, 0.315952], + [0.937990, 0.492667, 0.312575], + [0.940313, 0.497642, 0.309197], + [0.942598, 0.502639, 0.305816], + [0.944844, 0.507658, 0.302433], + [0.947051, 0.512699, 0.299049], + [0.949217, 0.517763, 0.295662], + [0.951344, 0.522850, 0.292275], + [0.953428, 0.527960, 0.288883], + [0.955470, 0.533093, 0.285490], + [0.957469, 0.538250, 0.282096], + [0.959424, 0.543431, 0.278701], + [0.961336, 0.548636, 0.275305], + [0.963203, 0.553865, 0.271909], + [0.965024, 0.559118, 0.268513], + [0.966798, 0.564396, 0.265118], + [0.968526, 0.569700, 0.261721], + [0.970205, 0.575028, 0.258325], + [0.971835, 0.580382, 0.254931], + [0.973416, 0.585761, 0.251540], + [0.974947, 0.591165, 0.248151], + [0.976428, 0.596595, 0.244767], + [0.977856, 0.602051, 0.241387], + [0.979233, 0.607532, 0.238013], + [0.980556, 0.613039, 0.234646], + [0.981826, 0.618572, 0.231287], + [0.983041, 0.624131, 0.227937], + [0.984199, 0.629718, 0.224595], + [0.985301, 0.635330, 0.221265], + [0.986345, 0.640969, 0.217948], + [0.987332, 0.646633, 0.214648], + [0.988260, 0.652325, 0.211364], + [0.989128, 0.658043, 0.208100], + [0.989935, 0.663787, 0.204859], + [0.990681, 0.669558, 0.201642], + [0.991365, 0.675355, 0.198453], + [0.991985, 0.681179, 0.195295], + [0.992541, 0.687030, 0.192170], + [0.993032, 0.692907, 0.189084], + [0.993456, 0.698810, 0.186041], + [0.993814, 0.704741, 0.183043], + [0.994103, 0.710698, 0.180097], + [0.994324, 0.716681, 0.177208], + [0.994474, 0.722691, 0.174381], + [0.994553, 0.728728, 0.171622], + [0.994561, 0.734791, 0.168938], + [0.994495, 0.740880, 0.166335], + [0.994355, 0.746995, 0.163821], + [0.994141, 0.753137, 0.161404], + [0.993851, 0.759304, 0.159092], + [0.993482, 0.765499, 0.156891], + [0.993033, 0.771720, 0.154808], + [0.992505, 0.777967, 0.152855], + [0.991897, 0.784239, 0.151042], + [0.991209, 0.790537, 0.149377], + [0.990439, 0.796859, 0.147870], + [0.989587, 0.803205, 0.146529], + [0.988648, 0.809579, 0.145357], + [0.987621, 0.815978, 0.144363], + [0.986509, 0.822401, 0.143557], + [0.985314, 0.828846, 0.142945], + [0.984031, 0.835315, 0.142528], + [0.982653, 0.841812, 0.142303], + [0.981190, 0.848329, 0.142279], + [0.979644, 0.854866, 0.142453], + [0.977995, 0.861432, 0.142808], + [0.976265, 0.868016, 0.143351], + [0.974443, 0.874622, 0.144061], + [0.972530, 0.881250, 0.144923], + [0.970533, 0.887896, 0.145919], + [0.968443, 0.894564, 0.147014], + [0.966271, 0.901249, 0.148180], + [0.964021, 0.907950, 0.149370], + [0.961681, 0.914672, 0.150520], + [0.959276, 0.921407, 0.151566], + [0.956808, 0.928152, 0.152409], + [0.954287, 0.934908, 0.152921], + [0.951726, 0.941671, 0.152925], + [0.949151, 0.948435, 0.152178], + [0.946602, 0.955190, 0.150328], + [0.944152, 0.961916, 0.146861], + [0.941896, 0.968590, 0.140956], + [0.940015, 0.975158, 0.131326]] + +_viridis_data = [[0.267004, 0.004874, 0.329415], + [0.268510, 0.009605, 0.335427], + [0.269944, 0.014625, 0.341379], + [0.271305, 0.019942, 0.347269], + [0.272594, 0.025563, 0.353093], + [0.273809, 0.031497, 0.358853], + [0.274952, 0.037752, 0.364543], + [0.276022, 0.044167, 0.370164], + [0.277018, 0.050344, 0.375715], + [0.277941, 0.056324, 0.381191], + [0.278791, 0.062145, 0.386592], + [0.279566, 0.067836, 0.391917], + [0.280267, 0.073417, 0.397163], + [0.280894, 0.078907, 0.402329], + [0.281446, 0.084320, 0.407414], + [0.281924, 0.089666, 0.412415], + [0.282327, 0.094955, 0.417331], + [0.282656, 0.100196, 0.422160], + [0.282910, 0.105393, 0.426902], + [0.283091, 0.110553, 0.431554], + [0.283197, 0.115680, 0.436115], + [0.283229, 0.120777, 0.440584], + [0.283187, 0.125848, 0.444960], + [0.283072, 0.130895, 0.449241], + [0.282884, 0.135920, 0.453427], + [0.282623, 0.140926, 0.457517], + [0.282290, 0.145912, 0.461510], + [0.281887, 0.150881, 0.465405], + [0.281412, 0.155834, 0.469201], + [0.280868, 0.160771, 0.472899], + [0.280255, 0.165693, 0.476498], + [0.279574, 0.170599, 0.479997], + [0.278826, 0.175490, 0.483397], + [0.278012, 0.180367, 0.486697], + [0.277134, 0.185228, 0.489898], + [0.276194, 0.190074, 0.493001], + [0.275191, 0.194905, 0.496005], + [0.274128, 0.199721, 0.498911], + [0.273006, 0.204520, 0.501721], + [0.271828, 0.209303, 0.504434], + [0.270595, 0.214069, 0.507052], + [0.269308, 0.218818, 0.509577], + [0.267968, 0.223549, 0.512008], + [0.266580, 0.228262, 0.514349], + [0.265145, 0.232956, 0.516599], + [0.263663, 0.237631, 0.518762], + [0.262138, 0.242286, 0.520837], + [0.260571, 0.246922, 0.522828], + [0.258965, 0.251537, 0.524736], + [0.257322, 0.256130, 0.526563], + [0.255645, 0.260703, 0.528312], + [0.253935, 0.265254, 0.529983], + [0.252194, 0.269783, 0.531579], + [0.250425, 0.274290, 0.533103], + [0.248629, 0.278775, 0.534556], + [0.246811, 0.283237, 0.535941], + [0.244972, 0.287675, 0.537260], + [0.243113, 0.292092, 0.538516], + [0.241237, 0.296485, 0.539709], + [0.239346, 0.300855, 0.540844], + [0.237441, 0.305202, 0.541921], + [0.235526, 0.309527, 0.542944], + [0.233603, 0.313828, 0.543914], + [0.231674, 0.318106, 0.544834], + [0.229739, 0.322361, 0.545706], + [0.227802, 0.326594, 0.546532], + [0.225863, 0.330805, 0.547314], + [0.223925, 0.334994, 0.548053], + [0.221989, 0.339161, 0.548752], + [0.220057, 0.343307, 0.549413], + [0.218130, 0.347432, 0.550038], + [0.216210, 0.351535, 0.550627], + [0.214298, 0.355619, 0.551184], + [0.212395, 0.359683, 0.551710], + [0.210503, 0.363727, 0.552206], + [0.208623, 0.367752, 0.552675], + [0.206756, 0.371758, 0.553117], + [0.204903, 0.375746, 0.553533], + [0.203063, 0.379716, 0.553925], + [0.201239, 0.383670, 0.554294], + [0.199430, 0.387607, 0.554642], + [0.197636, 0.391528, 0.554969], + [0.195860, 0.395433, 0.555276], + [0.194100, 0.399323, 0.555565], + [0.192357, 0.403199, 0.555836], + [0.190631, 0.407061, 0.556089], + [0.188923, 0.410910, 0.556326], + [0.187231, 0.414746, 0.556547], + [0.185556, 0.418570, 0.556753], + [0.183898, 0.422383, 0.556944], + [0.182256, 0.426184, 0.557120], + [0.180629, 0.429975, 0.557282], + [0.179019, 0.433756, 0.557430], + [0.177423, 0.437527, 0.557565], + [0.175841, 0.441290, 0.557685], + [0.174274, 0.445044, 0.557792], + [0.172719, 0.448791, 0.557885], + [0.171176, 0.452530, 0.557965], + [0.169646, 0.456262, 0.558030], + [0.168126, 0.459988, 0.558082], + [0.166617, 0.463708, 0.558119], + [0.165117, 0.467423, 0.558141], + [0.163625, 0.471133, 0.558148], + [0.162142, 0.474838, 0.558140], + [0.160665, 0.478540, 0.558115], + [0.159194, 0.482237, 0.558073], + [0.157729, 0.485932, 0.558013], + [0.156270, 0.489624, 0.557936], + [0.154815, 0.493313, 0.557840], + [0.153364, 0.497000, 0.557724], + [0.151918, 0.500685, 0.557587], + [0.150476, 0.504369, 0.557430], + [0.149039, 0.508051, 0.557250], + [0.147607, 0.511733, 0.557049], + [0.146180, 0.515413, 0.556823], + [0.144759, 0.519093, 0.556572], + [0.143343, 0.522773, 0.556295], + [0.141935, 0.526453, 0.555991], + [0.140536, 0.530132, 0.555659], + [0.139147, 0.533812, 0.555298], + [0.137770, 0.537492, 0.554906], + [0.136408, 0.541173, 0.554483], + [0.135066, 0.544853, 0.554029], + [0.133743, 0.548535, 0.553541], + [0.132444, 0.552216, 0.553018], + [0.131172, 0.555899, 0.552459], + [0.129933, 0.559582, 0.551864], + [0.128729, 0.563265, 0.551229], + [0.127568, 0.566949, 0.550556], + [0.126453, 0.570633, 0.549841], + [0.125394, 0.574318, 0.549086], + [0.124395, 0.578002, 0.548287], + [0.123463, 0.581687, 0.547445], + [0.122606, 0.585371, 0.546557], + [0.121831, 0.589055, 0.545623], + [0.121148, 0.592739, 0.544641], + [0.120565, 0.596422, 0.543611], + [0.120092, 0.600104, 0.542530], + [0.119738, 0.603785, 0.541400], + [0.119512, 0.607464, 0.540218], + [0.119423, 0.611141, 0.538982], + [0.119483, 0.614817, 0.537692], + [0.119699, 0.618490, 0.536347], + [0.120081, 0.622161, 0.534946], + [0.120638, 0.625828, 0.533488], + [0.121380, 0.629492, 0.531973], + [0.122312, 0.633153, 0.530398], + [0.123444, 0.636809, 0.528763], + [0.124780, 0.640461, 0.527068], + [0.126326, 0.644107, 0.525311], + [0.128087, 0.647749, 0.523491], + [0.130067, 0.651384, 0.521608], + [0.132268, 0.655014, 0.519661], + [0.134692, 0.658636, 0.517649], + [0.137339, 0.662252, 0.515571], + [0.140210, 0.665859, 0.513427], + [0.143303, 0.669459, 0.511215], + [0.146616, 0.673050, 0.508936], + [0.150148, 0.676631, 0.506589], + [0.153894, 0.680203, 0.504172], + [0.157851, 0.683765, 0.501686], + [0.162016, 0.687316, 0.499129], + [0.166383, 0.690856, 0.496502], + [0.170948, 0.694384, 0.493803], + [0.175707, 0.697900, 0.491033], + [0.180653, 0.701402, 0.488189], + [0.185783, 0.704891, 0.485273], + [0.191090, 0.708366, 0.482284], + [0.196571, 0.711827, 0.479221], + [0.202219, 0.715272, 0.476084], + [0.208030, 0.718701, 0.472873], + [0.214000, 0.722114, 0.469588], + [0.220124, 0.725509, 0.466226], + [0.226397, 0.728888, 0.462789], + [0.232815, 0.732247, 0.459277], + [0.239374, 0.735588, 0.455688], + [0.246070, 0.738910, 0.452024], + [0.252899, 0.742211, 0.448284], + [0.259857, 0.745492, 0.444467], + [0.266941, 0.748751, 0.440573], + [0.274149, 0.751988, 0.436601], + [0.281477, 0.755203, 0.432552], + [0.288921, 0.758394, 0.428426], + [0.296479, 0.761561, 0.424223], + [0.304148, 0.764704, 0.419943], + [0.311925, 0.767822, 0.415586], + [0.319809, 0.770914, 0.411152], + [0.327796, 0.773980, 0.406640], + [0.335885, 0.777018, 0.402049], + [0.344074, 0.780029, 0.397381], + [0.352360, 0.783011, 0.392636], + [0.360741, 0.785964, 0.387814], + [0.369214, 0.788888, 0.382914], + [0.377779, 0.791781, 0.377939], + [0.386433, 0.794644, 0.372886], + [0.395174, 0.797475, 0.367757], + [0.404001, 0.800275, 0.362552], + [0.412913, 0.803041, 0.357269], + [0.421908, 0.805774, 0.351910], + [0.430983, 0.808473, 0.346476], + [0.440137, 0.811138, 0.340967], + [0.449368, 0.813768, 0.335384], + [0.458674, 0.816363, 0.329727], + [0.468053, 0.818921, 0.323998], + [0.477504, 0.821444, 0.318195], + [0.487026, 0.823929, 0.312321], + [0.496615, 0.826376, 0.306377], + [0.506271, 0.828786, 0.300362], + [0.515992, 0.831158, 0.294279], + [0.525776, 0.833491, 0.288127], + [0.535621, 0.835785, 0.281908], + [0.545524, 0.838039, 0.275626], + [0.555484, 0.840254, 0.269281], + [0.565498, 0.842430, 0.262877], + [0.575563, 0.844566, 0.256415], + [0.585678, 0.846661, 0.249897], + [0.595839, 0.848717, 0.243329], + [0.606045, 0.850733, 0.236712], + [0.616293, 0.852709, 0.230052], + [0.626579, 0.854645, 0.223353], + [0.636902, 0.856542, 0.216620], + [0.647257, 0.858400, 0.209861], + [0.657642, 0.860219, 0.203082], + [0.668054, 0.861999, 0.196293], + [0.678489, 0.863742, 0.189503], + [0.688944, 0.865448, 0.182725], + [0.699415, 0.867117, 0.175971], + [0.709898, 0.868751, 0.169257], + [0.720391, 0.870350, 0.162603], + [0.730889, 0.871916, 0.156029], + [0.741388, 0.873449, 0.149561], + [0.751884, 0.874951, 0.143228], + [0.762373, 0.876424, 0.137064], + [0.772852, 0.877868, 0.131109], + [0.783315, 0.879285, 0.125405], + [0.793760, 0.880678, 0.120005], + [0.804182, 0.882046, 0.114965], + [0.814576, 0.883393, 0.110347], + [0.824940, 0.884720, 0.106217], + [0.835270, 0.886029, 0.102646], + [0.845561, 0.887322, 0.099702], + [0.855810, 0.888601, 0.097452], + [0.866013, 0.889868, 0.095953], + [0.876168, 0.891125, 0.095250], + [0.886271, 0.892374, 0.095374], + [0.896320, 0.893616, 0.096335], + [0.906311, 0.894855, 0.098125], + [0.916242, 0.896091, 0.100717], + [0.926106, 0.897330, 0.104071], + [0.935904, 0.898570, 0.108131], + [0.945636, 0.899815, 0.112838], + [0.955300, 0.901065, 0.118128], + [0.964894, 0.902323, 0.123941], + [0.974417, 0.903590, 0.130215], + [0.983868, 0.904867, 0.136897], + [0.993248, 0.906157, 0.143936]] + +from matplotlib.colors import ListedColormap + +cmaps = {} +for (name, data) in (('magma', _magma_data), + ('inferno', _inferno_data), + ('plasma', _plasma_data), + ('viridis', _viridis_data)): + + cmaps[name] = ListedColormap(data, name=name) + +magma = cmaps['magma'] +inferno = cmaps['inferno'] +plasma = cmaps['plasma'] +viridis = cmaps['viridis'] diff --git a/ice2ocn_fwflux.py b/ice2ocn_fwflux.py index d3ac4b0..cd66690 100644 --- a/ice2ocn_fwflux.py +++ b/ice2ocn_fwflux.py @@ -16,6 +16,9 @@ def ice2ocn_fwflux (file_path, tstep, fig_name): # Read data id = Dataset(file_path, 'r') + # I think there is a problem with salinity weighting here. + # There is a difference between kg/m^2/s of freshwater, and kg/m^2/s of + # salt. data_tmp = id.variables['fresh_ai'][tstep-1,:-15,:] - id.variables['fsalt_ai'][tstep-1,:-15,:]/1000.*60.*60.*24.*100. lon_tmp = id.variables['TLON'][:-15,:] lat_tmp = id.variables['TLAT'][:-15,:] diff --git a/ice_drift_seasonal.py b/ice_drift_seasonal.py new file mode 100644 index 0000000..d071718 --- /dev/null +++ b/ice_drift_seasonal.py @@ -0,0 +1,230 @@ +from numpy import * +from netCDF4 import Dataset, num2date +from matplotlib.pyplot import * +from rotate_vector_cice import * + +def ice_drift_seasonal (cice_file, save=False, fig_name=None): + + start_month = [2, 5, 8, 11] + end_month = [4, 7, 10, 1] + start_day = [1, 1, 1, 1] + end_day = [30, 31, 31, 31] + ndays_season = [89, 92, 92, 92] + season_names = ['FMA', 'MMJ', 'ASO', 'NDJ'] + figure_order = [1, 2, 4, 3] + r = 6.371e6 + deg2rad = pi/180.0 + block = 15 + + id = Dataset(cice_file, 'r') + lon_tmp = id.variables['TLON'][:-15,:] + lat_tmp = id.variables['TLAT'][:-15,:] + angle_tmp = id.variables['ANGLET'][:-15,:] + lon = ma.empty([size(lon_tmp,0), size(lon_tmp,1)+1]) + lat = ma.empty([size(lat_tmp,0), size(lat_tmp,1)+1]) + angle = ma.empty([size(angle_tmp,0), size(angle_tmp,1)+1]) + lon[:,:-1] = lon_tmp + lon[:,-1] = lon_tmp[:,0] + lat[:,:-1] = lat_tmp + lat[:,-1] = lat_tmp[:,0] + angle[:,:-1] = angle_tmp + angle[:,-1] = angle_tmp[:,0] + time_id = id.variables['time'] + time = num2date(time_id[:], units=time_id.units, calendar=time_id.calendar.lower()) + + end_t = -1 + for t in range(size(time)-1, -1, -1): + if time[t].month == start_month[0] and time[t].day in range(start_day[0], start_day[0]+5): + end_t = t + break + if end_t == -1: + print 'Error: ' + cice_file + ' does not contain a complete Feb-Jan period' + return + + start_t = -1 + for t in range(end_t-60, -1, -1): + if time[t].month == start_month[0] and time[t].day in range(start_day[0]+1, start_day[0]+6): + start_t = t + break + if start_t == -1: + print 'Error: ' + cice_file + ' does not contain a complete Feb-Jan period' + return + + leap_year = False + if mod(time[start_t].year, 4) == 0: + leap_year = True + if mod(time[start_t].year, 100) == 0: + leap_year = False + if mod(time[start_t].year, 400) == 0: + leap_year = True + if leap_year: + ndays_season[0] += 1 + + aice_tmp = ma.empty([4, size(lon_tmp,0), size(lon_tmp,1)]) + aice_tmp[:,:,:] = 0.0 + uxy_tmp = ma.empty([4, size(lon_tmp,0), size(lon_tmp,1)]) + uxy_tmp[:,:,:] = 0.0 + vxy_tmp = ma.empty([4, size(lon_tmp,0), size(lon_tmp,1)]) + vxy_tmp[:,:,:] = 0.0 + + for season in range(4): + season_days = 0 + next_season = mod(season+1, 4) + + start_t_season = -1 + for t in range(start_t, end_t+1): + if time[t].month == start_month[season] and time[t].day in range(start_day[season]+1, start_day[season]+6): + start_t_season = t + break + if start_t_season == -1: + print 'Error: could not find starting timestep for season ' + season_names[season] + return + + end_t_season = -1 + for t in range(start_t_season+1, end_t+1): + if time[t].month == start_month[next_season] and time[t].day in range(start_day[next_season], start_day[next_season]+5): + end_t_season = t + break + if end_t_season == -1: + print 'Error: could not find ending timestep for season ' + season_names[season] + return + + if time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season] + 5: + start_days = 5 + elif time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season] + 4: + start_days = 4 + elif time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season]+ 3: + start_days = 3 + elif time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season] + 2: + start_days = 2 + elif time[start_t_season].month == start_month[season] and time[start_t_season].day == start_day[season] + 1: + start_days = 1 + else: + print 'Error for season ' + season_names[season] + ': starting index is month ' + str(time[start_t_season].month) + ', day ' + str(time[start_t_season].day) + return + + aice_tmp[season,:,:] += id.variables['aice'][start_t_season,:-15,:]*start_days + uxy_tmp[season,:,:] += id.variables['uvel'][start_t_season,:-15,:]*start_days + vxy_tmp[season,:,:] += id.variables['vvel'][start_t_season,:-15,:]*start_days + season_days += start_days + + for t in range(start_t_season+1, end_t_season): + aice_tmp[season,:,:] += id.variables['aice'][t,:-15,:]*5 + uxy_tmp[season,:,:] += id.variables['uvel'][t,:-15,:]*5 + vxy_tmp[season,:,:] += id.variables['vvel'][t,:-15,:]*5 + season_days += 5 + + if time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season] + 4: + end_days = 1 + elif time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season] + 3: + end_days = 2 + elif time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season] + 2: + end_days = 3 + elif time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season] + 1: + end_days = 4 + elif time[end_t_season].month == start_month[next_season] and time[end_t_season].day == start_day[next_season]: + end_days = 5 + else: + print 'Error for season ' + season_names[season] + ': ending index is month ' + str(time[end_t_season].month) + ', day ' + str(time[end_t_season].day) + return + + aice_tmp[season,:,:] += id.variables['aice'][end_t_season,:-15,:]*end_days + uxy_tmp[season,:,:] += id.variables['uvel'][end_t_season,:-15,:]*end_days + vxy_tmp[season,:,:] += id.variables['vvel'][end_t_season,:-15,:]*end_days + + season_days += end_days + + if season_days != ndays_season[season]: + print 'Error: found ' + str(season_days) + ' days instead of ' + str(ndays_season[season]) + return + + aice_tmp[season,:,:] /= season_days + uxy_tmp[season,:,:] /= season_days + vxy_tmp[season,:,:] /= season_days + + id.close() + + aice = ma.empty([size(aice_tmp,0), size(aice_tmp,1), size(aice_tmp,2)+1]) + aice[:,:,:-1] = aice_tmp + aice[:,:,-1] = aice_tmp[:,:,0] + u_xy = ma.empty([size(uxy_tmp,0), size(uxy_tmp,1), size(uxy_tmp,2)+1]) + u_xy[:,:,:-1] = uxy_tmp + u_xy[:,:,-1] = uxy_tmp[:,:,0] + v_xy = ma.empty([size(vxy_tmp,0), size(vxy_tmp,1), size(vxy_tmp,2)+1]) + v_xy[:,:,:-1] = vxy_tmp + v_xy[:,:,-1] = vxy_tmp[:,:,0] + + u, v = rotate_vector_cice(u_xy, v_xy, angle) + speed = sqrt(u**2 + v**2) + theta = arctan2(v, u) + theta_circ = theta - lon*deg2rad + u_circ = speed*cos(theta_circ) + v_circ = speed*sin(theta_circ) + + x = -(lat+90)*cos(lon*deg2rad+pi/2) + y = (lat+90)*sin(lon*deg2rad+pi/2) + + size0 = int(ceil(size(x,0)/float(block))) + size1 = int(ceil((size(x,1)-1)/float(block))) + x_block = ma.empty([size0, size1]) + y_block = ma.empty([size0, size1]) + u_circ_block = ma.empty([4, size0, size1]) + v_circ_block = ma.empty([4, size0, size1]) + for season in range(4): + posn0 = range(0, size(x,0), block) + posn0.append(size(x,0)) + posn1 = range(0, size(x,1), block) + posn1.append(size(x,1)) + for j in range(size0): + for i in range(size1): + start0 = posn0[j] + end0 = posn0[j+1] + start1 = posn1[i] + end1 = posn1[i+1] + if season == 0: + x_block[j,i] = mean(x[start0:end0, start1:end1]) + y_block[j,i] = mean(y[start0:end0, start1:end1]) + u_circ_block[season,j,i] = mean(u_circ[season, start0:end0, start1:end1]) + v_circ_block[season,j,i] = mean(v_circ[season, start0:end0, start1:end1]) + + lev = linspace(0, 1, num=50) + bdry1 = -35 + bdry2 = 35 + bdry3 = -33 + bdry4 = 37 + + fig = figure(figsize=(16,12)) + for season in range(4): + ax = fig.add_subplot(2, 2, figure_order[season], aspect='equal') + img = contourf(x, y, aice[season,:,:], lev, cmap='jet') + q = quiver(x_block, y_block, u_circ_block[season,:,:], v_circ_block[season,:,:], color='black') + xlim([bdry1, bdry2]) + ylim([bdry3, bdry4]) + axis('off') + title(season_names[season], fontsize=24) + if season == 0: + cbaxes = fig.add_axes([0.07, 0.6, 0.02, 0.3]) + cbar = colorbar(img, ticks=arange(0,1+0.25,0.25), cax=cbaxes) + cbar.ax.tick_params(labelsize=16) + if season == 3: + quiverkey(q, 0.07, 0.3, 0.2, '20 cm/s', coordinates='figure', fontproperties={'size':16}) + suptitle('Sea ice concentration (1) and velocity (m/s)', fontsize=30) + subplots_adjust(wspace=0.025,hspace=0.1) + + if save: + fig.savefig(fig_name) + else: + fig.show() + + +if __name__ == "__main__": + + cice_file = raw_input("Path to CICE file, containing at least one complete Feb-Jan period: ") + action = raw_input("Save figure (s) or display on screen (d)? ") + if action == 's': + save = True + fig_name = raw_input("File name for figure: ") + elif action == 'd': + save = False + fig_name = None + ice_drift_seasonal(cice_file, save, fig_name) diff --git a/nsidc_aice_monthly.py b/nsidc_aice_monthly.py index 5866a77..8e6cbbd 100644 --- a/nsidc_aice_monthly.py +++ b/nsidc_aice_monthly.py @@ -41,7 +41,7 @@ def nsidc_aice_monthly (cice_file, month, save=False, fig_name=None): cice_lat[:,-1] = cice_lat_tmp[:,0] time_id = id.variables['time'] # Get the year, month, and day (all 1-based) for each output step - # These are 5-day averages marked with the last day's date. + # These are 5-day averages marked with the next day's date. cice_time = num2date(time_id[:], units=time_id.units, calendar=time_id.calendar.lower()) # Loop backwards through time indices to find the last one we care about @@ -49,11 +49,7 @@ def nsidc_aice_monthly (cice_file, month, save=False, fig_name=None): end_t = -1 # Missing value flag for t in range(size(cice_time)-1, -1, -1): # Note that cice_time[t].month is converted to 0-indexed - if cice_time[t].month-1 == month and cice_time[t].day == end_day[month]: - end_t = t - break - next_month = mod(month+1, 12) - if cice_time[t].month-1 == next_month and cice_time[t].day in range(start_day[next_month], start_day[next_month]+4): + if cice_time[t].month-1 == next_month and cice_time[t].day in range(start_day[next_month], start_day[next_month]+5): end_t = t break # Make sure we actually found it @@ -64,7 +60,7 @@ def nsidc_aice_monthly (cice_file, month, save=False, fig_name=None): # Continue looping backwards to find the first time index we care about start_t = -1 # Missing value flag for t in range(end_t, -1, -1): - if cice_time[t].month-1 == month and cice_time[t].day in range(start_day[month], start_day[month]+5): + if cice_time[t].month-1 == month and cice_time[t].day in range(start_day[month]+1, start_day[month]+6): start_t = t break # Make sure we actually found it @@ -74,14 +70,17 @@ def nsidc_aice_monthly (cice_file, month, save=False, fig_name=None): # Check if this is a leap year leap_year = False - if mod(cice_time[end_t].year, 4) == 0: + cice_year = cice_time[end_t].year + if month == 11: + cice_year = cice_time[start_t].year + if mod(cice_year, 4) == 0: # Years divisible by 4 are leap years leap_year = True - if mod(cice_time[end_t].year, 100) == 0: + if mod(cice_year, 100) == 0: # Unless they're also divisible by 100, in which case they aren't # leap years leap_year = False - if mod(cice_time[end_t].year, 400) == 0: + if mod(cice_year, 400) == 0: # Unless they're also divisible by 400, in which case they are # leap years after all leap_year = True @@ -96,19 +95,19 @@ def nsidc_aice_monthly (cice_file, month, save=False, fig_name=None): # Figure out how many of the 5 days averaged in start_t are actually within # this month - if cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month] + 4: + if cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month] + 5: # Starting day is in position 1 of 5; we care about all of them start_days = 5 - elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month] + 3: + elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month] + 4: # Starting day is in position 2 of 5; we care about the last 4 start_days = 4 - elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month]+ 2: + elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month]+ 3: # Starting day is in position 3 of 5; we care about the last 3 start_days = 3 - elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month] + 1: + elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month] + 2: # Starting day is in position 4 of 5; we care about the last 2 start_days = 2 - elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month]: + elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month] + 1: # Starting day is in position 5 of 5; we care about the last 1 start_days = 1 else: @@ -127,19 +126,19 @@ def nsidc_aice_monthly (cice_file, month, save=False, fig_name=None): # Figure out how many of the 5 days averaged in end_t are actually within # this month next_month = mod(month+1, 12) - if cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month] + 3: + if cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month] + 4: # Ending day is in position 1 of 5; we care about the first 1 end_days = 1 - elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month] + 2: + elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month] + 3: # Ending day is in position 2 of 5; we care about the first 2 end_days = 2 - elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month] + 1: + elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month] + 2: # Ending day is in position 3 of 5; we care about the first 3 end_days = 3 - elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month]: + elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month] + 1: # Ending day is in position 4 of 5; we care about the first 4 end_days = 4 - elif cice_time[end_t].month-1 == month and cice_time[end_t].day == end_day[month]: + elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month]: # Ending day is in position 5 of 5; we care about all 5 end_days = 5 else: diff --git a/nsidc_aice_seasonal.py b/nsidc_aice_seasonal.py index 006ac42..9ca0acd 100644 --- a/nsidc_aice_seasonal.py +++ b/nsidc_aice_seasonal.py @@ -47,17 +47,14 @@ def nsidc_aice_seasonal (cice_file, save=False, fig_name=None): cice_lat[:,-1] = cice_lat_tmp[:,0] time_id = id.variables['time'] # Get the year, month, and day (all 1-based) for each output step - # These are 5-day averages marked with the last day's date. + # These are 5-day averages marked with the next day's date. cice_time = num2date(time_id[:], units=time_id.units, calendar=time_id.calendar.lower()) # Loop backwards through time indices to find the last one we care about # (which contains 30 November in its averaging period) end_t = -1 # Missing value flag for t in range(size(cice_time)-1, -1, -1): - if cice_time[t].month == end_month[-1] and cice_time[t].day == end_day[-1]: - end_t = t - break - if cice_time[t].month == start_month[0] and cice_time[t].day in range(start_day[0], start_day[0]+4): + if cice_time[t].month == start_month[0] and cice_time[t].day in range(start_day[0], start_day[0]+5): end_t = t break # Make sure we actually found it @@ -69,7 +66,7 @@ def nsidc_aice_seasonal (cice_file, save=False, fig_name=None): # (which contains 1 December the previous year in its averaging period) start_t = -1 # Missing value flag for t in range(end_t-60, -1, -1): - if cice_time[t].month == start_month[0] and cice_time[t].day in range(start_day[0], start_day[0]+5): + if cice_time[t].month == start_month[0] and cice_time[t].day in range(start_day[0]+1, start_day[0]+6): start_t = t break # Make sure we actually found it @@ -77,7 +74,7 @@ def nsidc_aice_seasonal (cice_file, save=False, fig_name=None): print 'Error: ' + cice_file + ' does not contain a complete Dec-Nov period' return - # Check if end_t occurs on a leap year + # Check for leap years leap_year = False if mod(cice_time[end_t].year, 4) == 0: # Years divisible by 4 are leap years @@ -106,7 +103,7 @@ def nsidc_aice_seasonal (cice_file, save=False, fig_name=None): # Find starting timestep start_t_season = -1 for t in range(start_t, end_t+1): - if cice_time[t].month == start_month[season] and cice_time[t].day in range(start_day[season], start_day[season]+5): + if cice_time[t].month == start_month[season] and cice_time[t].day in range(start_day[season]+1, start_day[season]+6): start_t_season = t break # Make sure we actually found it @@ -117,10 +114,7 @@ def nsidc_aice_seasonal (cice_file, save=False, fig_name=None): # Find ending timestep end_t_season = -1 for t in range(start_t_season+1, end_t+1): - if cice_time[t].month == end_month[season] and cice_time[t].day == end_day[season]: - end_t_season = t - break - if cice_time[t].month == start_month[next_season] and cice_time[t].day in range(start_day[next_season], start_day[next_season]+4): + if cice_time[t].month == start_month[next_season] and cice_time[t].day in range(start_day[next_season]+1, start_day[next_season]+5): end_t_season = t break # Make sure we actually found it @@ -130,19 +124,19 @@ def nsidc_aice_seasonal (cice_file, save=False, fig_name=None): # Figure out how many of the 5 days averaged in start_t_season are # actually within this season - if cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 4: + if cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 5: # Starting day is in position 1 of 5; we care about all of them start_days = 5 - elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 3: + elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 4: # Starting day is in position 2 of 5; we care about the last 4 start_days = 4 - elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season]+ 2: + elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season]+ 3: # Starting day is in position 3 of 5; we care about the last 3 start_days = 3 - elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 1: + elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 2: # Starting day is in position 4 of 5; we care about the last 2 start_days = 2 - elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season]: + elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 1: # Starting day is in position 5 of 5; we care about the last 1 start_days = 1 else: @@ -160,19 +154,19 @@ def nsidc_aice_seasonal (cice_file, save=False, fig_name=None): # Figure out how many of the 5 days averaged in end_t_season are # actually within this season - if cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 3: + if cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 4: # Ending day is in position 1 of 5; we care about the first 1 end_days = 1 - elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 2: + elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 3: # Ending day is in position 2 of 5; we care about the first 2 end_days = 2 - elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 1: + elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 2: # Ending day is in position 3 of 5; we care about the first 3 end_days = 3 - elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season]: + elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 1: # Ending day is in position 4 of 5; we care about the first 4 end_days = 4 - elif cice_time[end_t_season].month == end_month[season] and cice_time[end_t_season].day == end_day[season]: + elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season]: # Ending day is in position 5 of 5; we care about all 5 end_days = 5 else: @@ -206,7 +200,7 @@ def nsidc_aice_seasonal (cice_file, save=False, fig_name=None): # Initialise seasonal averages of NSIDC data nsidc_data = ma.empty([4, size(nsidc_lon,0), size(nsidc_lon,1)]) - nsidc_data[:,:] = 0.0 + nsidc_data[:,:,:] = 0.0 # Process one season at a time for season in range(4): # Figure out which months we care about diff --git a/sose_seasonal.py b/sose_seasonal.py new file mode 100644 index 0000000..4343065 --- /dev/null +++ b/sose_seasonal.py @@ -0,0 +1,147 @@ +from netCDF4 import Dataset, num2date +from numpy import * +from matplotlib.pyplot import * +from calc_z import * + +def sose_seasonal (lon0, depth_bdry, save=False, fig_name=None): + + # Path to SOSE seasonal climatology file + sose_file = '../SOSE_seasonal_climatology.nc' + # Season names for titles + season_names = ['DJF', 'MAM', 'JJA', 'SON'] + + # Bounds on colour scale + temp_min = -2.5 + temp_max = 7.5 + temp_ticks = 2 + salt_min = 33.8 + salt_max = 34.8 + salt_ticks = 0.2 + + # Choose what to write on the title about longitude + if lon0 < 0: + lon_string = ' at ' + str(int(round(-lon0))) + r'$^{\circ}$W' + else: + lon_string = ' at ' + str(int(round(lon0))) + r'$^{\circ}$E' + # Edit longitude bounds to be from 0 to 360, to fit with ROMS convention + if lon0 < 0: + lon0 += 360 + + print 'Processing SOSE data' + + # Read grid and 3D data (already seasonally averaged) + id = Dataset(sose_file, 'r') + lon_sose = id.variables['longitude'][0,:] + lat_sose = id.variables['latitude'][:,0] + z_sose = id.variables['depth'][:] + temp_3d_sose = id.variables['temp'][:,:,:,:] + salt_3d_sose = id.variables['salt'][:,:,:,:] + + # Calculate zonal slices for each season + temp_sose = ma.empty([4, size(z_sose), size(lat_sose,0)]) + temp_sose[:,:,:] = 0.0 + salt_sose = ma.empty([4, size(z_sose), size(lat_sose,0)]) + salt_sose[:,:,:] = 0.0 + for season in range(4): + print 'Calculating zonal slices for ' + season_names[season] + temp_sose[season,:,:] = interp_lon_sose(temp_3d_sose[season,:,:,:], lon_sose, lon0) + salt_sose[season,:,:] = interp_lon_sose(salt_3d_sose[season,:,:,:], lon_sose, lon0) + + # Set colour levels + lev1 = linspace(temp_min, temp_max, num=50) + lev2 = linspace(salt_min, salt_max, num=50) + + # Choose southern boundary based on extent of SOSE grid + sbdry = -80 + nbdry = -30 + + # Plot + print 'Plotting' + fig = figure(figsize=(20,9)) + # Loop over seasons + for season in range(4): + # Temperature + fig.add_subplot(2, 4, season+1) + img = contourf(lat_sose, z_sose, temp_sose[season,:,:], lev1, cmap='jet', extend='both') + xlim([sbdry, nbdry]) + ylim([depth_bdry, 0]) + title('Temperature (' + season_names[season] + ')', fontsize=24) + if season == 0: + ylabel('depth (m)', fontsize=18) + if season == 3: + cbaxes1 = fig.add_axes([0.92, 0.55, 0.01, 0.3]) + cbar1 = colorbar(img, cax=cbaxes1, ticks=arange(temp_min, temp_max+temp_ticks, temp_ticks)) + cbar1.ax.tick_params(labelsize=16) + # Salinity + fig.add_subplot(2, 4, season+5) + img = contourf(lat_sose, z_sose, salt_sose[season,:,:], lev2, cmap='jet', extend='both') + xlim([sbdry, nbdry]) + ylim([depth_bdry, 0]) + title('Salinity (' + season_names[season] + ')', fontsize=24) + if season == 0: + ylabel('depth (m)', fontsize=18) + xlabel('Latitude', fontsize=18) + if season == 3: + cbaxes2 = fig.add_axes([0.92, 0.15, 0.01, 0.3]) + cbar2 = colorbar(img, cax=cbaxes2, ticks=arange(salt_min, salt_max+salt_ticks, salt_ticks)) + cbar2.ax.tick_params(labelsize=16) + # Add the main title + suptitle(lon_string, fontsize=30) + + # Finished + if save: + fig.savefig(fig_name) + else: + fig.show() + + +# Linearly interpolate SOSE data to the specified longitude. +# Input: +# data_3d = arary of data, dimension depth x lat x lon +# lon = 1D array of longitude values (between 0 and 360) +# lon0 = longitude to interpolate to (between 0 and 360) +# Output: +# data = array of data interpolated to lon0, dimension depth x lat +def interp_lon_sose (data_3d, lon, lon0): + + if lon0 < lon[0] or lon0 > lon[-1]: + # Special case: lon0 on periodic boundary + # Be careful with mod 360 here + iw = size(lon)-1 + ie = 0 + # Calculate difference between lon0 and lon[iw], mod 360 if necessary + dlon_num = lon0 - lon[iw] + if dlon_num < -300: + dlon_num += 360 + # Calculate difference between lon[ie] and lon[iw], mod 360 + dlon_den = lon[ie] - lon[iw] + 360 + else: + # General case + # Find the first index eastwards of lon0 + ie = nonzero(lon > lon0)[0][0] + # The index before it will be the last index westward of lon0 + iw = ie - 1 + dlon_num = lon0 - lon[iw] + dlon_den = lon[ie] - lon[iw] + # Coefficients for interpolation + coeff1 = dlon_num/dlon_den + coeff2 = 1 - coeff1 + + # Interpolate + data = coeff1*data_3d[:,:,ie] + coeff2*data_3d[:,:,iw] + return data + + +# Command-line interface +if __name__ == "__main__": + + lon0 = float(raw_input("Enter longitude (-180 to 180): ")) + depth_bdry = -1*float(raw_input("Deepest depth to plot (positive, metres): ")) + action = raw_input("Save figure (s) or display in window (d)? ") + if action == 's': + save = True + fig_name = raw_input("File name for figure: ") + elif action == 'd': + save = False + fig_name = None + sose_seasonal(lon0, depth_bdry, save, fig_name) diff --git a/ssflux_tamura_monthly.py b/ssflux_tamura_monthly.py new file mode 100644 index 0000000..d9440f3 --- /dev/null +++ b/ssflux_tamura_monthly.py @@ -0,0 +1,209 @@ +from numpy import * +from netCDF4 import Dataset, num2date +from matplotlib.pyplot import * + +def ssflux_tamura_monthly (cice_file, month, save=False, fig_name=None): + + tamura_head = '/short/m68/kaa561/tamura_fluxes/Tamura_ssflux_' + tamura_tail = '_monthly.nc' + start_day = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] + end_day = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] + month_name = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October', 'November', 'December'] + deg2rad = pi/180.0 + rho_fw = 1000.0 + rho_sw = 1026.0 + mps_to_cmpday = 8.64e6 + + id = Dataset(cice_file, 'r') + cice_lon_tmp = id.variables['TLON'][:-15,:] + cice_lat_tmp = id.variables['TLAT'][:-15,:] + cice_lon = ma.empty([size(cice_lon_tmp,0), size(cice_lon_tmp,1)+1]) + cice_lat = ma.empty([size(cice_lat_tmp,0), size(cice_lat_tmp,1)+1]) + cice_lon[:,:-1] = cice_lon_tmp + cice_lon[:,-1] = cice_lon_tmp[:,0] + cice_lat[:,:-1] = cice_lat_tmp + cice_lat[:,-1] = cice_lat_tmp[:,0] + time_id = id.variables['time'] + cice_time = num2date(time_id[:], units=time_id.units, calendar='standard') + + next_month = mod(month+1, 12) + prev_month = mod(month-1, 12) + + end_t = -1 + for t in range(size(cice_time)-1, -1, -1): + if cice_time[t].month-1 == next_month and cice_time[t].day in range(start_day[next_month], start_day[next_month]+5): + end_t = t + break + if end_t == -1: + print 'Error: ' + cice_file + ' does not contain a complete ' + month_name[month] + return + + start_t = -1 + for t in range(end_t, -1, -1): + if cice_time[t].month-1 == month and cice_time[t].day in range(start_day[month]+1, start_day[month]+6): + start_t = t + break + if start_t == -1: + print 'Error: ' + cice_file + ' does not contain a complete ' + month_name[month] + return + + leap_year = False + cice_year = cice_time[end_t].year + if month == 11: + cice_year = cice_time[start_t].year + if mod(cice_year, 4) == 0: + leap_year = True + if mod(cice_year, 100) == 0: + leap_year = False + if mod(cice_.year, 400) == 0: + leap_year = True + if leap_year: + end_day[1] = 29 + + cice_data_tmp = ma.empty(shape(cice_lon_tmp)) + cice_data_tmp[:,:] = 0.0 + num_days = 0 + + if cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month]+5: + start_days = 5 + elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month]+4: + start_days = 4 + elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month]+3: + start_days = 3 + elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month]+2: + start_days = 2 + elif cice_time[start_t].month-1 == month and cice_time[start_t].day == start_day[month]+1: + start_days = 1 + else: + print 'Error: starting index is month ' + str(cice_time[start_t].month) + ', day ' + str(cice_time[start_t].day) + return + + fresh_ai = id.variables['fresh_ai'][start_t,:-15,:] + fsalt_ai = id.variables['fsalt_ai'][start_t,:-15,:] + sss = id.variables['sss'][start_t,:-15,:] + cice_data_tmp += -1/rho_fw*(fresh_ai*rho_sw/mps_to_cmpday - fsalt_ai*1e3/sss)*start_days + num_days += start_days + + for t in range(start_t+1, end_t): + fresh_ai = id.variables['fresh_ai'][t,:-15,:] + fsalt_ai = id.variables['fsalt_ai'][t,:-15,:] + sss = id.variables['sss'][t,:-15,:] + cice_data_tmp += -1/rho_fw*(fresh_ai*rho_sw/mps_to_cmpday - fsalt_ai*1e3/sss)*5 + num_days += 5 + + if cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month]+4: + end_days = 1 + elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month]+3: + end_days = 2 + elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month]+2: + end_days = 3 + elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month]+1: + end_days = 4 + elif cice_time[end_t].month-1 == next_month and cice_time[end_t].day == start_day[next_month]: + end_days = 5 + else: + print 'Error: ending index is month ' + str(cice_time[end_t].month) + ', day ' + str(cice_time[end_t].day) + return + + fresh_ai = id.variables['fresh_ai'][end_t,:-15,:] + fsalt_ai = id.variables['fsalt_ai'][end_t,:-15,:] + sss = id.variables['sss'][end_t,:-15,:] + cice_data_tmp += -1/rho_fw*(fresh_ai*rho_sw/mps_to_cmpday - fsalt_ai*1e3/sss)*end_days + num_days += end_days + + if num_days != end_day[month]: + print 'Error: found ' + str(num_days) + ' days instead of ' + str(end_day[month]) + return + + id.close() + cice_data_tmp /= num_days + cice_data_tmp *= 1e6 + + cice_data = ma.empty([size(cice_data_tmp,0), size(cice_data_tmp,1)+1]) + cice_data[:,:-1] = cice_data_tmp + cice_data[:,-1] = cice_data_tmp[:,0] + + id = Dataset(tamura_head + str(cice_year) + tamura_tail, 'r') + tamura_lon = id.variables['longitude'][:,:] + tamura_lat = id.variables['latitude'][:,:] + tamura_data = id.variables['ssflux'][month,:,:] + id.close() + tamura_data = ma.masked_where(isnan(tamura_data), tamura_data) + tamura_data *= 1e6 + + cice_x = -(cice_lat+90)*cos(cice_lon*deg2rad+pi/2) + cice_y = (cice_lat+90)*sin(cice_lon*deg2rad+pi/2) + tamura_x = -(tamura_lat+90)*cos(tamura_lon*deg2rad+pi/2) + tamura_y = (tamura_lat+90)*sin(tamura_lon*deg2rad+pi/2) + + bound = 1.0 #max(amax(abs(cice_data)), amax(abs(tamura_data))) + lev = linspace(-bound, bound, num=50) + bdry1 = max(amin(cice_x), amin(tamura_x)) + bdry2 = min(amax(cice_x), amax(tamura_x)) + bdry3 = max(amin(cice_y), amin(tamura_y)) + bdry4 = min(amax(cice_y), amax(tamura_y)) + + fig = figure(figsize=(20,9)) + fig.add_subplot(1,2,1, aspect='equal') + contourf(tamura_x, tamura_y, tamura_data, lev, cmap='RdYlBu_r') + title('Tamura', fontsize=24) + xlim([bdry1, bdry2]) + ylim([bdry3, bdry4]) + axis('off') + fig.add_subplot(1,2,2, aspect='equal') + img = contourf(cice_x, cice_y, cice_data, lev, cmap='RdYlBu_r') + title('CICE', fontsize=24) + xlim([bdry1, bdry2]) + ylim([bdry3, bdry4]) + axis('off') + cbaxes = fig.add_axes([0.3, 0.04, 0.4, 0.04]) + cbar = colorbar(img, orientation='horizontal', ticks=arange(-1,1+0.5, 0.5), cax=cbaxes, extend='both') + cbar.ax.tick_params(labelsize=18) + suptitle(r'Ice-to-ocean salt flux (10$^{-6}$ kg/m$^2$/s, ' + month_name[month] + ' ' + str(cice_year), fontsize=30) + + if save: + fig.savefig(fig_name) + else: + fig.show() + + +if __name__ == "__main__": + + cice_file = raw_input("Path to CICE file: ") + month = int(raw_input("Month number (1-12): "))-1 + action = raw_input("Save figure (s) or display on screen (d)? ") + if action == 's': + save = True + fig_name = raw_input("File name for figure: ") + elif action == 'd': + save = False + fig_name = None + ssflux_tamura_monthly(cice_file, month, save, fig_name) + + while True: + repeat = raw_input("Make another plot (y/n)? ") + if repeat == 'y': + while True: + changes = raw_input("Enter a parameter to change: (1) file path, (2) month number, (3) save/display; or enter to continue: ") + if len(changes) == 0: + break + else: + if int(changes) == 1: + cice_file = raw_input("Path to CICE file: ") + elif int(changes) == 2: + month = int(raw_input("Month number (1-12): "))-1 + elif int(changes) == 3: + save = not save + if save: + fig_name = raw_input("File name for figure: ") + ssflux_tamura_monthly(cice_file, month, save, fig_name) + else: + break + + + + + + + + diff --git a/ssflux_tamura_seasonal.py b/ssflux_tamura_seasonal.py new file mode 100644 index 0000000..504e1b4 --- /dev/null +++ b/ssflux_tamura_seasonal.py @@ -0,0 +1,233 @@ +from numpy import * +from netCDF4 import Dataset, num2date +from matplotlib.pyplot import * + +def ssflux_tamura_seasonal (cice_file, save=False, fig_name=None): + + tamura_head = '/short/m68/kaa561/tamura_fluxes/Tamura_ssflux_' + tamura_tail = '_monthly.nc' + start_month = [12, 3, 6, 9] + end_month = [2, 5, 8, 11] + start_day = [1, 1, 1, 1] + end_day = [28, 31, 31, 30] + ndays_season = [90, 92, 92, 91] + ndays_month = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] + season_names = ['DJF', 'MAM', 'JJA', 'SON'] + deg2rad = pi/180.0 + rho_fw = 1000.0 + rho_sw = 1026.0 + mps_to_cmpday = 8.64e6 + + id = Dataset(cice_file, 'r') + cice_lon_tmp = id.variables['TLON'][:-15,:] + cice_lat_tmp = id.variables['TLAT'][:-15,:] + cice_lon = ma.empty([size(cice_lon_tmp,0), size(cice_lon_tmp,1)+1]) + cice_lat = ma.empty([size(cice_lat_tmp,0), size(cice_lat_tmp,1)+1]) + cice_lon[:,:-1] = cice_lon_tmp + cice_lon[:,-1] = cice_lon_tmp[:,0] + cice_lat[:,:-1] = cice_lat_tmp + cice_lat[:,-1] = cice_lat_tmp[:,0] + time_id = id.variables['time'] + cice_time = num2date(time_id[:], units=time_id.units, calendar=time_id.calendar.lower()) + + end_t = -1 + for t in range(size(cice_time)-1, -1, -1): + if cice_time[t].month == start_month[0] and cice_time[t].day in range(start_day[0], start_day[0]+5): + end_t = t + break + if end_t == -1: + print 'Error: ' + cice_file + ' does not contain a complete Dec-Nov period' + return + + start_t = -1 + for t in range(end_t-60, -1, -1): + if cice_time[t].month == start_month[0] and cice_time[t].day in range(start_day[0]+1, start_day[0]+6): + start_t = t + break + if start_t == -1: + print 'Error: ' + cice_file + ' does not contain a complete Dec-Nov period' + return + + leap_year = False + cice_year = cice_time[end_t].year + if mod(cice_year, 4) == 0: + leap_year = True + if mod(cice_year, 100) == 0: + leap_year = False + if mod(cice_year, 400) == 0: + leap_year = True + if leap_year: + end_day[0] += 1 + ndays_season[0] += 1 + + cice_data_tmp = ma.empty([4, size(cice_lon_tmp,0), size(cice_lon_tmp,1)]) + cice_data_tmp[:,:,:] = 0.0 + + for season in range(4): + season_days = 0 + next_season = mod(season+1, 4) + + start_t_season = -1 + for t in range(start_t, end_t+1): + if cice_time[t].month == start_month[season] and cice_time[t].day in range(start_day[season]+1, start_day[season]+6): + start_t_season = t + break + if start_t_season == -1: + print 'Error: could not find starting timestep for season ' + season_names[season] + return + + end_t_season = -1 + for t in range(start_t_season+1, end_t+1): + if cice_time[t].month == start_month[next_season] and cice_time[t].day in range(start_day[next_season], start_day[next_season]+5): + end_t_season = t + break + if end_t_season == -1: + print 'Error: could not find ending timestep for season ' + season_names[season] + return + + if cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 5: + start_days = 5 + elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 4: + start_days = 4 + elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season]+ 3: + start_days = 3 + elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 2: + start_days = 2 + elif cice_time[start_t_season].month == start_month[season] and cice_time[start_t_season].day == start_day[season] + 1: + start_days = 1 + else: + print 'Error for season ' + season_names[season] + ': starting index is month ' + str(cice_time[start_t_season].month) + ', day ' + str(cice_time[start_t_season].day) + return + + fresh_ai = id.variables['fresh_ai'][start_t_season,:-15,:] + rain_ai = id.variables['rain_ai'][start_t_season,:-15,:] + fsalt_ai = id.variables['fsalt_ai'][start_t_season,:-15,:] + sss = id.variables['sss'][start_t_season,:-15,:] + cice_data_tmp[season,:,:] += -1/rho_fw*((fresh_ai-rain_ai)*rho_sw/mps_to_cmpday - fsalt_ai*1e3/sss)*start_days + season_days += start_days + + for t in range(start_t_season+1, end_t_season): + fresh_ai = id.variables['fresh_ai'][t,:-15,:] + rain_ai = id.variables['rain_ai'][t,:-15,:] + fsalt_ai = id.variables['fsalt_ai'][t,:-15,:] + sss = id.variables['sss'][t,:-15,:] + cice_data_tmp[season,:,:] += -1/rho_fw*((fresh_ai-rain_ai)*rho_sw/mps_to_cmpday - fsalt_ai*1e3/sss)*5 + season_days += 5 + + if cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 4: + end_days = 1 + elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 3: + end_days = 2 + elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 2: + end_days = 3 + elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season] + 1: + end_days = 4 + elif cice_time[end_t_season].month == start_month[next_season] and cice_time[end_t_season].day == start_day[next_season]: + end_days = 5 + else: + print 'Error for season ' + season_names[season] + ': ending index is month ' + str(cice_time[end_t_season].month) + ', day ' + str(cice_time[end_t_season].day) + return + + fresh_ai = id.variables['fresh_ai'][end_t_season,:-15,:] + rain_ai = id.variables['rain_ai'][end_t_season,:-15,:] + fsalt_ai = id.variables['fsalt_ai'][end_t_season,:-15,:] + sss = id.variables['sss'][end_t_season,:-15,:] + cice_data_tmp[season,:,:] += -1/rho_fw*((fresh_ai-rain_ai)*rho_sw/mps_to_cmpday - fsalt_ai*1e3/sss)*end_days + season_days += end_days + + if season_days != ndays_season[season]: + print 'Error: found ' + str(season_days) + ' days instead of ' + str(ndays_season[season]) + return + + cice_data_tmp[season,:,:] /= season_days + + id.close() + cice_data_tmp *= 1e6 + + cice_data = ma.empty([size(cice_data_tmp,0), size(cice_data_tmp,1), size(cice_data_tmp,2)+1]) + cice_data[:,:,:-1] = cice_data_tmp + cice_data[:,:,-1] = cice_data_tmp[:,:,0] + + id = Dataset(tamura_head + str(cice_year) + tamura_tail, 'r') + tamura_lon = id.variables['longitude'][:,:] + tamura_lat = id.variables['latitude'][:,:] + id.close() + + tamura_data = ma.empty([4, size(tamura_lon,0), size(tamura_lon,1)]) + tamura_data[:,:,:] = 0.0 + for season in range(4): + if season == 0: + tamura_months = [11, 0, 1] + elif season == 1: + tamura_months = [2, 3, 4] + elif season == 2: + tamura_months = [5, 6, 7] + elif season == 3: + tamura_months = [8, 9, 10] + season_days = 0 + + for month in tamura_months: + if season == 0 and month == 11: + id = Dataset(tamura_head + str(cice_year-1) + tamura_tail, 'r') + else: + id = Dataset(tamura_head + str(cice_year) + tamura_tail, 'r') + tamura_data_tmp = id.variables['ssflux'][month,:,:] + id.close() + tamura_data_tmp = ma.masked_where(isnan(tamura_data_tmp), tamura_data_tmp) + tamura_data[season,:,:] += tamura_data_tmp*ndays_month[month] + season_days += ndays_month[month] + tamura_data[season,:,:] /= season_days + + tamura_data *= 1e6 + + cice_x = -(cice_lat+90)*cos(cice_lon*deg2rad+pi/2) + cice_y = (cice_lat+90)*sin(cice_lon*deg2rad+pi/2) + tamura_x = -(tamura_lat+90)*cos(tamura_lon*deg2rad+pi/2) + tamura_y = (tamura_lat+90)*sin(tamura_lon*deg2rad+pi/2) + + lev = linspace(-0.5, 0.5, num=50) + bdry1 = -35 #max(amin(cice_x), amin(tamura_x)) + bdry2 = 35 #min(amax(cice_x), amax(tamura_x)) + bdry3 = -33 #max(amin(cice_y), amin(tamura_y)) + bdry4 = 37 #min(amax(cice_y), amax(tamura_y)) + + fig = figure(figsize=(20,9)) + for season in range(4): + ax = fig.add_subplot(2, 4, season+1, aspect='equal') + contourf(tamura_x, tamura_y, tamura_data[season,:,:], lev, cmap='RdYlBu_r', extend='both') + if season == 0: + text(-39, 0, 'Tamura', fontsize=24, ha='right') + title(season_names[season], fontsize=24) + xlim([bdry1, bdry2]) + ylim([bdry3, bdry4]) + axis('off') + ax = fig.add_subplot(2, 4, season+5, aspect='equal') + img = contourf(cice_x, cice_y, cice_data[season,:,:], lev, cmap='RdYlBu_r', extend='both') + if season == 0: + text(-39, 0, 'CICE', fontsize=24, ha='right') + xlim([bdry1, bdry2]) + ylim([bdry3, bdry4]) + axis('off') + cbaxes = fig.add_axes([0.25, 0.04, 0.5, 0.02]) + cbar = colorbar(img, orientation='horizontal', ticks=arange(-0.5,0.5+0.25, 0.25), cax=cbaxes) + cbar.ax.tick_params(labelsize=16) + suptitle(r'Ice-to-ocean salt flux (10$^{-6}$ kg/m$^2$/s)', fontsize=30) + subplots_adjust(wspace=0.025,hspace=0.025) + + if save: + fig.savefig(fig_name) + else: + fig.show() + + +if __name__ == "__main__": + + cice_file = raw_input("Path to CICE file, containing at least one complete Dec-Nov period: ") + action = raw_input("Save figure (s) or display on screen (d)? ") + if action == 's': + save = True + fig_name = raw_input("File name for figure: ") + elif action == 'd': + save = False + fig_name = None + ssflux_tamura_seasonal(cice_file, save, fig_name) diff --git a/temp_salt_seasonal.py b/temp_salt_seasonal.py index d5f2cf5..f863a89 100644 --- a/temp_salt_seasonal.py +++ b/temp_salt_seasonal.py @@ -42,9 +42,9 @@ def temp_salt_seasonal (file_path, lon0, depth_bdry, save=False, fig_name=None): temp_min = -2.5 temp_max = 7.5 temp_ticks = 2 - salt_min = 34.3 + salt_min = 33.8 salt_max = 34.8 - salt_ticks = 0.1 + salt_ticks = 0.2 # Choose what to write on the title about longitude if lon0 < 0: diff --git a/uv_vectorplot.py b/uv_vectorplot.py index 5928a85..cf4ef24 100644 --- a/uv_vectorplot.py +++ b/uv_vectorplot.py @@ -54,33 +54,29 @@ def uv_vectorplot (grid_path, file_path, tstep, depth_key, save=False, fig_name= # Calculate speed for the background filled contour plot speed = sqrt(u_rho**2 + v_rho**2) - # Calculate X and Y coordinates for plotting circumpolar projection - X = -(lat+90)*cos(lon*deg2rad+pi/2) - Y = (lat+90)*sin(lon*deg2rad+pi/2) + # Calculate x and y coordinates for plotting circumpolar projection + x = -(lat+90)*cos(lon*deg2rad+pi/2) + y = (lat+90)*sin(lon*deg2rad+pi/2) - # Calculate velocity components in spherical coordinate space - # (just differentiate and rearrange spherical coordinate transformations) - dlon_dt = u_rho/(r*cos(lat*deg2rad)*deg2rad) - dlat_dt = v_rho/(r*deg2rad) - # Calculate velocity components in X-Y space (just differentiate and - # rearrange equations for X and Y above) - dX_dt = -dlat_dt*cos(lon*deg2rad+pi/2) + (lat+90)*sin(lon*deg2rad+pi/2)*dlon_dt*deg2rad - dY_dt = dlat_dt*sin(lon*deg2rad+pi/2) + (lat+90)*cos(lon*deg2rad+pi/2)*dlon_dt*deg2rad + theta = arctan2(v_rho, u_rho) + theta_circ = theta - lon*deg2rad + u_circ = speed*cos(theta_circ) + v_circ = speed*sin(theta_circ) # Average X, Y, dX_dt, and dY_dt over block x block intervals # Calculate number of blocks - size0 = int(ceil(size(X,0)/float(block))) - size1 = int(ceil((size(X,1)-1)/float(block))) + size0 = int(ceil(size(x,0)/float(block))) + size1 = int(ceil((size(x,1)-1)/float(block))) # Set up arrays for averaged fields - X_block = ma.empty([size0, size1]) - Y_block = ma.empty([size0, size1]) - dX_dt_block = ma.empty([size0, size1]) - dY_dt_block = ma.empty([size0, size1]) + x_block = ma.empty([size0, size1]) + y_block = ma.empty([size0, size1]) + u_circ_block = ma.empty([size0, size1]) + v_circ_block = ma.empty([size0, size1]) # Set up arrays containing boundary indices - posn0 = range(0, size(X,0), block) - posn0.append(size(X,0)) - posn1 = range(0, size(X,1), block) - posn1.append(size(X,1)) + posn0 = range(0, size(x,0), block) + posn0.append(size(x,0)) + posn1 = range(0, size(x,1), block) + posn1.append(size(x,1)) # Double loop to average each block (can't find a more efficient way to do # this) for j in range(size0): @@ -89,21 +85,21 @@ def uv_vectorplot (grid_path, file_path, tstep, depth_key, save=False, fig_name= end0 = posn0[j+1] start1 = posn1[i] end1 = posn1[i+1] - X_block[j,i] = mean(X[start0:end0, start1:end1]) - Y_block[j,i] = mean(Y[start0:end0, start1:end1]) - dX_dt_block[j,i] = mean(dX_dt[start0:end0, start1:end1]) - dY_dt_block[j,i] = mean(dY_dt[start0:end0, start1:end1]) + x_block[j,i] = mean(x[start0:end0, start1:end1]) + y_block[j,i] = mean(y[start0:end0, start1:end1]) + u_circ_block[j,i] = mean(u_circ[start0:end0, start1:end1]) + v_circ_block[j,i] = mean(v_circ[start0:end0, start1:end1]) # Make the plot fig = figure(figsize=(16,12)) fig.add_subplot(1,1,1, aspect='equal') # Contour speed values at every point # Use pastel colour map so overlaid vectors will show up - contourf(X, Y, speed, 50, cmap='Paired') + contourf(x, y, speed, 50, cmap='Paired') cbar = colorbar() cbar.ax.tick_params(labelsize=20) # Add vectors for each block - quiver(X_block, Y_block, dX_dt_block, dY_dt_block, color='black') + quiver(x_block, y_block, u_circ_block, v_circ_block, color='black') if depth_key == 1: title('Surface velocity (m/s)', fontsize=30) elif depth_key == 2: diff --git a/zonal_plot.py b/zonal_plot.py index 4d5abb5..ff878d5 100644 --- a/zonal_plot.py +++ b/zonal_plot.py @@ -172,7 +172,7 @@ def zonal_plot (file_path, var_name, tstep, lon_key, lon0, lon_bounds, depth_min # There is land everywhere at the northern boundary # Show the first 2 degrees of this land mask lat_max = max(lat[:,j_max]) + 2 -# lat_max = -65 +# lat_max = -75 xlim([lat_min, lat_max]) ylim([depth_min, 0])