From 302343ecbaa45cf6a63e4a4ba9e6345e653f2e43 Mon Sep 17 00:00:00 2001 From: Kaitlin Naughten Date: Wed, 9 Aug 2017 14:23:45 +1000 Subject: [PATCH] Tweaks to ice shelf figures --- file_guide.txt | 12 ++++ mip_aice_minmax_nsidc.py | 18 ++++-- mip_grid_res.py | 9 ++- mip_hi_seasonal.py | 21 ++++--- mip_iceshelf_figures.py | 83 +++++++++++++++------------ mip_mld.py | 35 +++++++++--- mip_seasonal_cycle.py | 117 +++++++++++++++++++++++++++++++++++++++ timeseries_dpt.py | 6 +- timeseries_massloss.py | 6 +- timeseries_seaice.py | 6 +- 10 files changed, 248 insertions(+), 65 deletions(-) create mode 100644 mip_seasonal_cycle.py diff --git a/file_guide.txt b/file_guide.txt index f64ea7f..9a16385 100644 --- a/file_guide.txt +++ b/file_guide.txt @@ -1436,6 +1436,18 @@ mip_iceshelf_figures.py: This is the giant monster script to generate 8 GridSpec arguments, colourbar axes, etc. so that the positioning looks okay. +mip_seasonal_cycle.py: Calculate the average amplitude of the seasonal cycle in + total basal mass loss over 2002-2016, in MetROMS, low-res + FESOM, and high-res FESOM. Print results to the screen. + To run: First make sure you have run + timeseries_massloss.py for the entire MetROMS + simulation, and the equivalent fesomtools + script for both FESOM simulations. Then open + python or ipython and type + "run mip_seasonal_cycle.py". The script will + prompt you for the paths to the + timeseries_massloss logfiles for each simulation. + ***UTILITY FUNCTIONS*** diff --git a/mip_aice_minmax_nsidc.py b/mip_aice_minmax_nsidc.py index c6f895e..53ee4d3 100644 --- a/mip_aice_minmax_nsidc.py +++ b/mip_aice_minmax_nsidc.py @@ -306,7 +306,8 @@ def mip_aice_minmax_nsidc (cice_file, cice_log, fesom_mesh_path, fesom_output_di img = pcolor(nsidc_x, nsidc_y, nsidc_feb, vmin=0, vmax=1, cmap='jet') xlim([bdry1, bdry2]) ylim([bdry3, bdry4]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) title('NSIDC', fontsize=24) text(-39, 0, 'February', fontsize=24, ha='right') # MetROMS, February @@ -314,7 +315,8 @@ def mip_aice_minmax_nsidc (cice_file, cice_log, fesom_mesh_path, fesom_output_di img = pcolor(cice_x, cice_y, cice_feb, vmin=0, vmax=1, cmap='jet') xlim([bdry1, bdry2]) ylim([bdry3, bdry4]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) title('MetROMS', fontsize=24) # FESOM, February ax = subplot(gs1[0, 2], aspect='equal') @@ -325,7 +327,8 @@ def mip_aice_minmax_nsidc (cice_file, cice_log, fesom_mesh_path, fesom_output_di ax.add_collection(img) xlim([bdry1, bdry2]) ylim([bdry3, bdry4]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) title('FESOM (high-res)', fontsize=24) # Main title text(-170, 47, 'a) Sea ice concentration ('+str(start_year)+'-'+str(end_year)+' average)', fontsize=30) @@ -334,14 +337,16 @@ def mip_aice_minmax_nsidc (cice_file, cice_log, fesom_mesh_path, fesom_output_di img = pcolor(nsidc_x, nsidc_y, nsidc_sep, vmin=0, vmax=1, cmap='jet') xlim([bdry1, bdry2]) ylim([bdry3, bdry4]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) text(-39, 0, 'September', fontsize=24, ha='right') # MetROMS, September ax = subplot(gs1[1, 1], aspect='equal') img = pcolor(cice_x, cice_y, cice_sep, vmin=0, vmax=1, cmap='jet') xlim([bdry1, bdry2]) ylim([bdry3, bdry4]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) # FESOM, September ax = subplot(gs1[1, 2], aspect='equal') img = PatchCollection(patches, cmap='jet') @@ -351,7 +356,8 @@ def mip_aice_minmax_nsidc (cice_file, cice_log, fesom_mesh_path, fesom_output_di ax.add_collection(img) xlim([bdry1, bdry2]) ylim([bdry3, bdry4]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) # Add a colourbar at the bottom cbaxes = fig.add_axes([0.26, 0.04, 0.3, 0.04]) cbar = colorbar(img, orientation='horizontal', ticks=arange(0,1+0.25,0.25), cax=cbaxes) diff --git a/mip_grid_res.py b/mip_grid_res.py index 5dfc7e2..1787905 100644 --- a/mip_grid_res.py +++ b/mip_grid_res.py @@ -69,7 +69,8 @@ def mip_grid_res (roms_grid_file, fesom_mesh_low, fesom_mesh_high, save=False, f pcolor(roms_x, roms_y, roms_res, vmin=limits[0], vmax=limits[1], cmap='jet') xlim([-lat_max, lat_max]) ylim([-lat_max, lat_max]) - axis('off') + ax1.set_xticks([]) + ax1.set_yticks([]) title('a) MetROMS', fontsize=28) # FESOM low-res ax2 = fig.add_subplot(1,3,2, aspect='equal') @@ -80,7 +81,8 @@ def mip_grid_res (roms_grid_file, fesom_mesh_low, fesom_mesh_high, save=False, f ax2.add_collection(img_low) xlim([-lat_max, lat_max]) ylim([-lat_max, lat_max]) - axis('off') + ax2.set_xticks([]) + ax2.set_yticks([]) title('b) FESOM low-res', fontsize=28) # FESOM high-res ax3 = fig.add_subplot(1,3,3, aspect='equal') @@ -91,7 +93,8 @@ def mip_grid_res (roms_grid_file, fesom_mesh_low, fesom_mesh_high, save=False, f ax3.add_collection(img_high) xlim([-lat_max, lat_max]) ylim([-lat_max, lat_max]) - axis('off') + ax3.set_xticks([]) + ax3.set_yticks([]) title('c) FESOM high-res', fontsize=28) cbaxes = fig.add_axes([0.92, 0.2, 0.01, 0.6]) cbar = colorbar(img_high, cax=cbaxes, extend='max', ticks=arange(limits[0], limits[1]+5, 5)) diff --git a/mip_hi_seasonal.py b/mip_hi_seasonal.py index f71e911..350deb7 100644 --- a/mip_hi_seasonal.py +++ b/mip_hi_seasonal.py @@ -18,8 +18,11 @@ # "fesomtools" repository def mip_hi_seasonal (cice_seasonal_file, fesom_mesh_path, fesom_seasonal_file): - # Northern boundary of plot 50S - nbdry = -50 + 90 + # Boundaries on plot (under polar coordinate transformation) + x_min = -36.25 + x_max = 36.25 + y_min = -34.5 + y_max = 38 # Degrees to radians conversion factor deg2rad = pi/180.0 # FESOM parameters @@ -89,9 +92,10 @@ def mip_hi_seasonal (cice_seasonal_file, fesom_mesh_path, fesom_seasonal_file): if season == 0: text(-43, 0, 'MetROMS', fontsize=24, ha='right') title(season_names[season], fontsize=24) - xlim([-nbdry, nbdry]) - ylim([-nbdry, nbdry]) - axis('off') + xlim([x_min, x_max]) + ylim([y_min, y_max]) + ax.set_xticks([]) + ax.set_yticks([]) # FESOM ax = fig.add_subplot(2, 4, season+5, aspect='equal') img = PatchCollection(patches, cmap='jet') @@ -99,9 +103,10 @@ def mip_hi_seasonal (cice_seasonal_file, fesom_mesh_path, fesom_seasonal_file): img.set_clim(vmin=bounds[0], vmax=bounds[1]) img.set_edgecolor('face') ax.add_collection(img) - xlim([-nbdry, nbdry]) - ylim([-nbdry, nbdry]) - axis('off') + xlim([x_min, x_max]) + ylim([y_min, y_max]) + ax.set_xticks([]) + ax.set_yticks([]) if season == 0: text(-43, 0, 'FESOM', fontsize=24, ha='right') text(-43, -10, '(high-res)', fontsize=24,ha='right') diff --git a/mip_iceshelf_figures.py b/mip_iceshelf_figures.py index 1c74213..0b67386 100644 --- a/mip_iceshelf_figures.py +++ b/mip_iceshelf_figures.py @@ -499,7 +499,7 @@ gs_a.update(left=0.05, right=0.9, bottom=0.735, top=0.89, wspace=0.05) cbaxes_tmp = fig.add_axes([0.91, 0.755, 0.025, 0.12]) cbar_ticks = arange(0, 6+3, 3) -plot_melt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_a, cbaxes_tmp, cbar_ticks, [0.5, 3, 4.5], 'a', 1.25, [-60, -40], [-80]) +plot_melt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_a, cbaxes_tmp, cbar_ticks, [0.5, 3, 4.5], 'a', 1.25, [-80, -60, -40], [-80]) # Velocity x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr = make_vectors(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, 20, 20) gs_b = GridSpec(1,3) @@ -511,7 +511,7 @@ gs_c = GridSpec(1,3) gs_c.update(left=0.05, right=0.9, bottom=0.375, top=0.53, wspace=0.05) cbaxes_tmp = fig.add_axes([0.91, 0.395, 0.025, 0.12]) -cbar_ticks = arange(-2.6, -1.8+0.4, 0.4) +cbar_ticks = arange(-2.6, -1.0+0.8, 0.8) plot_bwtemp(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_c, cbaxes_tmp, cbar_ticks, 'c') # Bottom water salinity gs_d = GridSpec(1,3) @@ -541,7 +541,7 @@ gs_a.update(left=0.11, right=0.9, bottom=0.74, top=0.89, wspace=0.05) cbaxes_tmp = fig.add_axes([0.91, 0.77, 0.025, 0.1]) cbar_ticks = arange(0, 4+2, 2) -plot_melt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_a, cbaxes_tmp, cbar_ticks, [1, 2, 3], 'a', 1.35, [0, 30], [-70]) +plot_melt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_a, cbaxes_tmp, cbar_ticks, [1, 2, 3], 'a', 1.35, [-10, 20], [-70]) # Velocity x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr = make_vectors(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, 40, 20) gs_b = GridSpec(1,3) @@ -630,13 +630,13 @@ gs_b = GridSpec(1,3) gs_b.update(left=0.05, right=0.9, bottom=0.435, top=0.655, wspace=0.05) cbaxes_tmp = fig.add_axes([0.91, 0.475, 0.025, 0.15]) -cbar_ticks = arange(-1.8, -0.8+0.5, 0.5) +cbar_ticks = arange(-1.5, 0.5+1, 1) plot_bwtemp(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_b, cbaxes_tmp, cbar_ticks, 'b') # Bottom water salinity gs_c = GridSpec(1,3) gs_c.update(left=0.05, right=0.9, bottom=0.19, top=0.41, wspace=0.05) cbaxes_tmp = fig.add_axes([0.91, 0.23, 0.025, 0.15]) -cbar_ticks = arange(34.2, 34.5+0.1, 0.1) +cbar_ticks = arange(34.2, 34.6+0.2, 0.2) plot_bwsalt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_c, cbaxes_tmp, cbar_ticks, 'c') # Velocity for Totten x_min_tmp = 20 @@ -665,14 +665,14 @@ gs_a.update(left=0.1, right=0.9, bottom=0.735, top=0.89, wspace=0.05) cbaxes_tmp = fig.add_axes([0.91, 0.765, 0.025, 0.1]) cbar_ticks = arange(0, 8+4, 4) -plot_melt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_a, cbaxes_tmp, cbar_ticks, [0.5, 2, 4], 'a', 1.25, [180, -160], [-80]) +plot_melt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_a, cbaxes_tmp, cbar_ticks, [0.5, 2, 4], 'a', 1.25, [180, -140], [-80]) # Velocity x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr = make_vectors(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, 20, 15) gs_b = GridSpec(1,3) gs_b.update(left=0.1, right=0.9, bottom=0.56, top=0.715, wspace=0.05) cbaxes_tmp = fig.add_axes([0.91, 0.59, 0.025, 0.1]) cbar_ticks = arange(0, 0.18+0.09, 0.09) -plot_velavg(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_b, cbaxes_tmp, cbar_ticks, x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr, 'b', arrow_scale=0.5, arrow_headwidth=8, arrow_headlength=9) +plot_velavg(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_b, cbaxes_tmp, cbar_ticks, x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr, 'b', arrow_scale=0.5, arrow_headwidth=7, arrow_headlength=8) # Draft gs_c = GridSpec(1,3) gs_c.update(left=0.1, right=0.9, bottom=0.385, top=0.54, wspace=0.05) @@ -696,7 +696,7 @@ gs_e.update(left=0.1, right=0.9, bottom=0.05, top=0.18, wspace=0.05) cbaxes_tmp2 = fig.add_axes([0.91, 0.065, 0.025, 0.1]) cbar_ticks2 = arange(34.4, 35+0.3, 0.3) -plot_zonal_ts(180, lat_min, lat_max, lat_ticks, lat_labels, depth_min, depth_max, depth_ticks, depth_labels, gs_d, gs_e, cbaxes_tmp1, cbaxes_tmp2, cbar_ticks1, cbar_ticks2, 'c', 'd') +plot_zonal_ts(180, lat_min, lat_max, lat_ticks, lat_labels, depth_min, depth_max, depth_ticks, depth_labels, gs_d, gs_e, cbaxes_tmp1, cbaxes_tmp2, cbar_ticks1, cbar_ticks2, 'd', 'e') suptitle('Ross Sea', fontsize=30) fig.show() fig.savefig('ross.png') @@ -718,13 +718,13 @@ gs_b = GridSpec(1,3) gs_b.update(left=0.05, right=0.9, bottom=0.435, top=0.655, wspace=0.05) cbaxes_tmp = fig.add_axes([0.91, 0.475, 0.025, 0.15]) -cbar_ticks = arange(-1.6, 0+0.8, 0.8) +cbar_ticks = arange(-1.6, 0.8+0.8, 0.8) plot_bwtemp(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_b, cbaxes_tmp, cbar_ticks, 'b') # Bottom water salinity gs_c = GridSpec(1,3) gs_c.update(left=0.05, right=0.9, bottom=0.19, top=0.41, wspace=0.05) cbaxes_tmp = fig.add_axes([0.91, 0.23, 0.025, 0.15]) -cbar_ticks = arange(34.1, 34.5+0.2, 0.2) +cbar_ticks = arange(34.1, 34.7+0.2, 0.2) plot_bwsalt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_c, cbaxes_tmp, cbar_ticks, 'c') # Velocity for PIG x_min_tmp = -15.6 @@ -758,13 +758,13 @@ gs_b = GridSpec(1,3) gs_b.update(left=0.05, right=0.9, bottom=0.435, top=0.655, wspace=0.05) cbaxes_tmp = fig.add_axes([0.91, 0.475, 0.025, 0.15]) -cbar_ticks = arange(-2, 0+1, 1) +cbar_ticks = arange(-1.5, 0.5+1, 1) plot_bwtemp(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_b, cbaxes_tmp, cbar_ticks, 'b') # Bottom water salinity gs_c = GridSpec(1,3) gs_c.update(left=0.05, right=0.9, bottom=0.19, top=0.41, wspace=0.05) cbaxes_tmp = fig.add_axes([0.91, 0.23, 0.025, 0.15]) -cbar_ticks = arange(33.6, 34.4+0.4, 0.4) +cbar_ticks = arange(33.8, 34.6+0.4, 0.4) plot_bwsalt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_c, cbaxes_tmp, cbar_ticks, 'c') # Velocity for George VI x_min_tmp = -18.75 @@ -791,21 +791,21 @@ # Melt gs_a = GridSpec(1,3) gs_a.update(left=0.05, right=0.9, bottom=0.61, top=0.84, wspace=0.05) -cbaxes_tmp = fig.add_axes([0.91, 0.7, 0.025, 0.15]) +cbaxes_tmp = fig.add_axes([0.91, 0.65, 0.025, 0.15]) cbar_ticks = arange(0, 8+4, 4) plot_melt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_a, cbaxes_tmp, cbar_ticks, [0.5, 2, 4], 'a', 1.3, [-60], [-70]) # Velocity x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr = make_vectors(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, 18, 18) gs_b = GridSpec(1,3) gs_b.update(left=0.05, right=0.9, bottom=0.33, top=0.56, wspace=0.05) -cbaxes_tmp = fig.add_axes([0.91, 0.4, 0.025, 0.15]) +cbaxes_tmp = fig.add_axes([0.91, 0.37, 0.025, 0.15]) cbar_ticks = arange(0, 0.12+0.06, 0.06) -plot_velavg(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_b, cbaxes_tmp, cbar_ticks, x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr, 'b', loc_string='George VI Ice Shelf', arrow_scale=0.7, arrow_headwidth=8, arrow_headlength=9) +plot_velavg(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_b, cbaxes_tmp, cbar_ticks, x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr, 'b', arrow_scale=0.7, arrow_headwidth=8, arrow_headlength=9) # Bottom water temperature gs_c = GridSpec(1,3) gs_c.update(left=0.05, right=0.9, bottom=0.05, top=0.28, wspace=0.05) -cbaxes_tmp = fig.add_axes([0.91, 0.1, 0.025, 0.15]) -cbar_ticks = arange(-1.8, -1+0.4, 0.4) +cbaxes_tmp = fig.add_axes([0.91, 0.09, 0.025, 0.15]) +cbar_ticks = arange(-1.8, 0+0.9, 0.9) plot_bwtemp(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_c, cbaxes_tmp, cbar_ticks, 'c') suptitle('Larsen Ice Shelves', fontsize=30) fig.show() @@ -987,7 +987,8 @@ def plot_draft (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter): pcolor(roms_x, roms_y, roms_draft, vmin=var_min, vmax=var_max, cmap='jet') xlim([x_min, x_max]) ylim([y_min, y_max]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) # FESOM low-res ax = subplot(gs[0,1], aspect='equal') @@ -1005,7 +1006,8 @@ def plot_draft (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter): ax.add_collection(overlay) xlim([x_min, x_max]) ylim([y_min, y_max]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) # Main title title(letter + ') Ice shelf draft (m)', fontsize=20) @@ -1022,7 +1024,8 @@ def plot_draft (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter): ax.add_collection(overlay) xlim([x_min, x_max]) ylim([y_min, y_max]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) # Colourbar on the right cbar = colorbar(img, cax=cbaxes, ticks=cbar_ticks) @@ -1093,12 +1096,13 @@ def plot_melt (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, change_points # Now shade the data pcolor(roms_x, roms_y, roms_melt, vmin=var_min, vmax=var_max, cmap=mf_cmap) # Overlay longitudes - contour(roms_x, roms_y, roms_lon, lon_lines, colors='black', linestyles='dashed') + contour(roms_x, roms_y, roms_lon, lon_lines, colors='black', linestyles='dotted') # Overlay latitudes - contour(roms_x, roms_y, roms_lat, lat_lines, colors='black', linestyles='dashed') + contour(roms_x, roms_y, roms_lat, lat_lines, colors='black', linestyles='dotted') xlim([x_min, x_max]) ylim([y_min, y_max]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) # Melt rate is always at the top, so add model labels text(0.5, y0, 'MetROMS', fontsize=18, horizontalalignment='center', transform=ax.transAxes) @@ -1118,7 +1122,8 @@ def plot_melt (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, change_points ax.add_collection(overlay) xlim([x_min, x_max]) ylim([y_min, y_max]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) title(letter + ') Ice shelf melt rate (m/y)', fontsize=20) text(0.5, y0, 'FESOM (low-res)', fontsize=18, horizontalalignment='center', transform=ax.transAxes) @@ -1135,7 +1140,8 @@ def plot_melt (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, change_points ax.add_collection(overlay) xlim([x_min, x_max]) ylim([y_min, y_max]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) # Colourbar on the right cbar = colorbar(img, cax=cbaxes, ticks=cbar_ticks) text(0.5, y0, 'FESOM (high-res)', fontsize=18, horizontalalignment='center', transform=ax.transAxes) @@ -1173,7 +1179,8 @@ def plot_bwtemp (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter): contour(roms_x, roms_y, zice_contour, levels=[min_zice], colors=('black')) xlim([x_min, x_max]) ylim([y_min, y_max]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) # FESOM low-res ax = subplot(gs[0,1], aspect='equal') @@ -1194,7 +1201,8 @@ def plot_bwtemp (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter): ax.add_collection(fesom_contours_lr) xlim([x_min, x_max]) ylim([y_min, y_max]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) # Main title title(letter + r') Bottom water temperature ($^{\circ}$C)', fontsize=20) @@ -1213,7 +1221,8 @@ def plot_bwtemp (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter): ax.add_collection(fesom_contours_hr) xlim([x_min, x_max]) ylim([y_min, y_max]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) # Colourbar on the right cbar = colorbar(img, cax=cbaxes, ticks=cbar_ticks) @@ -1250,7 +1259,8 @@ def plot_bwsalt (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter): contour(roms_x, roms_y, zice_contour, levels=[min_zice], colors=('black')) xlim([x_min, x_max]) ylim([y_min, y_max]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) # FESOM low-res ax = subplot(gs[0,1], aspect='equal') @@ -1271,7 +1281,8 @@ def plot_bwsalt (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter): ax.add_collection(fesom_contours_lr) xlim([x_min, x_max]) ylim([y_min, y_max]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) # Main title title(letter + ') Bottom water salinity (psu)', fontsize=20) @@ -1290,7 +1301,8 @@ def plot_bwsalt (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter): ax.add_collection(fesom_contours_hr) xlim([x_min, x_max]) ylim([y_min, y_max]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) # Colourbar on the right cbar = colorbar(img, cax=cbaxes, ticks=cbar_ticks) @@ -1334,7 +1346,8 @@ def plot_velavg (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, x_centres, quiver(x_centres, y_centres, roms_ubin, roms_vbin, scale=arrow_scale, headwidth=arrow_headwidth, headlength=arrow_headlength, color='black') xlim([x_min, x_max]) ylim([y_min, y_max]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) # FESOM low-res ax = subplot(gs[0,1], aspect='equal') @@ -1354,7 +1367,8 @@ def plot_velavg (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, x_centres, quiver(x_centres, y_centres, fesom_ubin_lr, fesom_vbin_lr, scale=arrow_scale, headwidth=arrow_headwidth, headlength=arrow_headlength, color='black') xlim([x_min, x_max]) ylim([y_min, y_max]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) # Main title if loc_string is None: title(letter + ') Vertically averaged ocean velocity (m/s)', fontsize=20) @@ -1375,7 +1389,8 @@ def plot_velavg (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, x_centres, quiver(x_centres, y_centres, fesom_ubin_hr, fesom_vbin_hr, scale=arrow_scale, headwidth=arrow_headwidth, headlength=arrow_headlength, color='black') xlim([x_min, x_max]) ylim([y_min, y_max]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) # Colourbar on the right cbar = colorbar(img, cax=cbaxes, ticks=cbar_ticks) diff --git a/mip_mld.py b/mip_mld.py index 79f79ff..4219ef0 100644 --- a/mip_mld.py +++ b/mip_mld.py @@ -48,6 +48,11 @@ def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file # Maximum for colour scale in each season max_bound_summer = 150 max_bound_winter = 600 + # Longitude labels for first panel + lon_ticks = array([-120, -60, 60, 120]) + lat_ticks = array([-28, -25, -25, -28]) + lon_labels = [r'120$^{\circ}$W', r'60$^{\circ}$W', r'60$^{\circ}$E', r'120$^{\circ}$E'] + lon_rot = [-60, 60, -60, 60] print 'Processing MetROMS:' print 'Reading grid' @@ -60,6 +65,9 @@ def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file # Polar coordinates for plotting roms_x = -(roms_lat+90)*cos(roms_lon*deg2rad+pi/2) roms_y = (roms_lat+90)*sin(roms_lon*deg2rad+pi/2) + # Longitude labels + x_ticks = -(lat_ticks+90)*cos(lon_ticks*deg2rad+pi/2) + y_ticks = (lat_ticks+90)*sin(lon_ticks*deg2rad+pi/2) # Get a 3D array of z-coordinates; sc_r and Cs_r are unused in this script roms_z, sc_r, Cs_r = calc_z(roms_h, roms_zice, theta_s, theta_b, hc, N) # Make depth positive @@ -150,28 +158,35 @@ def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file title(season_names[0], fontsize=24) xlim([-nbdry1, nbdry1]) ylim([-nbdry1, nbdry1]) - axis('off') + # Add longitude labels + for i in range(size(x_ticks)): + text(x_ticks[i], y_ticks[i], lon_labels[i], ha='center', rotation=lon_rot[i], fontsize=14) + ax.set_xticks([]) + ax.set_yticks([]) # MetROMS, summer, continental shelf ax = fig.add_subplot(2, 4, 2, aspect='equal') pcolor(roms_x, roms_y, roms_mld[0,:,:], vmin=0, vmax=max_bound_summer, cmap='jet') title(season_names[0], fontsize=24) xlim([-nbdry2, nbdry2]) ylim([-nbdry2, nbdry2]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) # MetROMS, winter, ACC ax = fig.add_subplot(2, 4, 3, aspect='equal') pcolor(roms_x, roms_y, roms_mld[2,:,:], vmin=0, vmax=max_bound_winter, cmap='jet') title(season_names[2], fontsize=24) xlim([-nbdry1, nbdry1]) ylim([-nbdry1, nbdry1]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) # MetROMS, winter, continental shelf ax = fig.add_subplot(2, 4, 4, aspect='equal') pcolor(roms_x, roms_y, roms_mld[2,:,:], vmin=0, vmax=max_bound_winter, cmap='jet') title(season_names[2], fontsize=24) xlim([-nbdry2, nbdry2]) ylim([-nbdry2, nbdry2]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) # FESOM, summer, ACC ax = fig.add_subplot(2, 4, 5, aspect='equal') img = PatchCollection(patches, cmap='jet') @@ -181,7 +196,8 @@ def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file ax.add_collection(img) xlim([-nbdry1, nbdry1]) ylim([-nbdry1, nbdry1]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) text(-62, 0, 'FESOM', fontsize=24, ha='right') text(-62, -10, '(high-res)', fontsize=24, ha='right') # FESOM, summer, continental shelf @@ -193,7 +209,8 @@ def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file ax.add_collection(img) xlim([-nbdry2, nbdry2]) ylim([-nbdry2, nbdry2]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) # Add a horizontal colorbar for summer cbaxes = fig.add_axes([0.2, 0.04, 0.2, 0.02]) cbar = colorbar(img, orientation='horizontal', cax=cbaxes, extend='max', ticks=arange(0, max_bound_summer+50, 50)) @@ -207,7 +224,8 @@ def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file ax.add_collection(img) xlim([-nbdry1, nbdry1]) ylim([-nbdry1, nbdry1]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) # FESOM, winter, continental shelf ax = fig.add_subplot(2, 4, 8, aspect='equal') img = PatchCollection(patches, cmap='jet') @@ -217,7 +235,8 @@ def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file ax.add_collection(img) xlim([-nbdry2, nbdry2]) ylim([-nbdry2, nbdry2]) - axis('off') + ax.set_xticks([]) + ax.set_yticks([]) # Add a horizontal colorbar for winter cbaxes = fig.add_axes([0.6, 0.04, 0.2, 0.02]) cbar = colorbar(img, orientation='horizontal', cax=cbaxes, extend='max', ticks=arange(0, max_bound_winter+200, 200)) diff --git a/mip_seasonal_cycle.py b/mip_seasonal_cycle.py new file mode 100644 index 0000000..5e62cdc --- /dev/null +++ b/mip_seasonal_cycle.py @@ -0,0 +1,117 @@ +from numpy import * + +# Calculate the average amplitude of the seasonal cycle in total basal mass +# loss over 2002-2016, in MetROMS, low-res FESOM, and high-res FESOM. Print +# results to the screen. +def mip_seasonal_cycle (roms_logfile, fesom_logfile_lr, fesom_logfile_hr): + + # Year simulations start + year_start = 1992 + # Years to avearge over + calc_start = 2002 + calc_end = 2016 + # Number of output steps per year in FESOM + peryear = 365/5 + + # Read ROMS logfile + roms_time = [] + f = open(roms_logfile, 'r') + # Skip the first line (header for time array) + f.readline() + for line in f: + try: + roms_time.append(float(line)) + except(ValueError): + # Reached the header for the next variable + break + # Read total mass loss + roms_massloss = [] + for line in f: + try: + roms_massloss.append(float(line)) + except(ValueError): + # Reached the header for the first individual ice shelf + break + f.close() + # Add start year to ROMS time array + roms_time = array(roms_time) + year_start + # Calculate amplitude for each year + roms_amplitude = [] + for year in range(calc_start, calc_end): + # Find the first index after the beginning of this year + t_start = nonzero(roms_time >= year)[0][0] + # Find the first index after the beginning of next year + tmp2 = nonzero(roms_time >= year+1)[0] + if len(tmp2)==0: + # No such index, but we might not have run out of data + # eg simulation that ends on 31 December 2016 + t_end = len(roms_time) + else: + t_end = tmp2[0] + # Get min and max between these bounds + curr_min = amin(roms_massloss[t_start:t_end]) + curr_max = amax(roms_massloss[t_start:t_end]) + # Save amplitude + roms_amplitude.append(curr_max-curr_min) + # Calculate average amplitude + roms_val = mean(array(roms_amplitude)) + print 'Average amplitude of seasonal cycle in total basal mass loss: ' + print 'ROMS: ' + str(roms_val) + + # Read FESOM logfiles + # Low-res + f = open(fesom_logfile_lr, 'r') + # Skip the first line (header) + f.readline() + # Read total mass loss + fesom_massloss_lr = [] + for line in f: + try: + fesom_massloss_lr.append(float(line)) + except(ValueError): + # Reached the header for the first individual ice shelf + break + f.close() + # Calculate amplitude for each year + fesom_amplitude_lr = [] + for year in range(calc_start, calc_end): + t_start = peryear*(year-year_start) + t_end = peryear*(year+1-year_start) + # Get min and max between these bounds + curr_min = amin(fesom_massloss_lr[t_start:t_end]) + curr_max = amax(fesom_massloss_lr[t_start:t_end]) + # Save amplitude + fesom_amplitude_lr.append(curr_max-curr_min) + # Calculate average amplitude + fesom_val_lr = mean(array(fesom_amplitude_lr)) + print 'FESOM low-res: ' + str(fesom_val_lr) + + # High-res + f = open(fesom_logfile_hr, 'r') + f.readline() + fesom_massloss_hr = [] + for line in f: + try: + fesom_massloss_hr.append(float(line)) + except(ValueError): + break + f.close() + fesom_amplitude_hr = [] + for year in range(calc_start, calc_end): + t_start = peryear*(year-year_start) + t_end = peryear*(year+1-year_start) + curr_min = amin(fesom_massloss_hr[t_start:t_end]) + curr_max = amax(fesom_massloss_hr[t_start:t_end]) + fesom_amplitude_hr.append(curr_max-curr_min) + fesom_val_hr = mean(array(fesom_amplitude_hr)) + print 'FESOM high-res: ' + str(fesom_val_hr) + + +# Command-line interface +if __name__ == "__main__": + + roms_logfile = raw_input("Path to ROMS logfile from timeseries_massloss.py: ") + fesom_logfile_lr = raw_input("Path to FESOM low-res logfile from timeseries_massloss.py: ") + fesom_logfile_hr = raw_input("Path to FESOM high-res logfile from timeseries_massloss.py: ") + mip_seasonal_cycle(roms_logfile, fesom_logfile_lr, fesom_logfile_hr) + diff --git a/timeseries_dpt.py b/timeseries_dpt.py index 12768ab..0644ed1 100644 --- a/timeseries_dpt.py +++ b/timeseries_dpt.py @@ -13,7 +13,9 @@ # log_path = path to log file (if it exists, previously calculated values will # be read from it; regardless, it will be overwritten with all # calculated values following computation) -def timeseries_dpt (grid_path, file_path, log_path): +# add_years = optional number of years to add to time array (multiple of 14 +# for repeating 1992-2005 spinup) +def timeseries_dpt (grid_path, file_path, log_path, add_years=0): # Radius of the Earth in metres r = 6.371e6 @@ -56,7 +58,7 @@ def timeseries_dpt (grid_path, file_path, log_path): print 'Reading data' id = Dataset(file_path, 'r') # Read time values and convert from seconds to years - new_time = id.variables['ocean_time'][:]/(60*60*24*365.25) + new_time = id.variables['ocean_time'][:]/(60*60*24*365.25) + add_years num_time = size(new_time) # Concatenate with time values from log file for t in range(num_time): diff --git a/timeseries_massloss.py b/timeseries_massloss.py index c7addc0..6c95c74 100644 --- a/timeseries_massloss.py +++ b/timeseries_massloss.py @@ -12,7 +12,9 @@ # log_path = path to log file (if it exists, previously calculated values will # be read from it; regardless, it will be overwritten with all # calculated values following computation) -def timeseries_massloss (file_path, log_path): +# add_years = optional number of years to add to time array (multiple of 14 +# for repeating 1992-2005 spinup) +def timeseries_massloss (file_path, log_path, add_years=0): # Titles and figure names for each ice shelf names = ['All Ice Shelves', 'Larsen D Ice Shelf', 'Larsen C Ice Shelf', 'Wilkins & George VI & Stange Ice Shelves', 'Ronne-Filchner Ice Shelf', 'Abbot Ice Shelf', 'Pine Island Glacier Ice Shelf', 'Thwaites Ice Shelf', 'Dotson Ice Shelf', 'Getz Ice Shelf', 'Nickerson Ice Shelf', 'Sulzberger Ice Shelf', 'Mertz Ice Shelf', 'Totten & Moscow University Ice Shelves', 'Shackleton Ice Shelf', 'West Ice Shelf', 'Amery Ice Shelf', 'Prince Harald Ice Shelf', 'Baudouin & Borchgrevink Ice Shelves', 'Lazarev Ice Shelf', 'Nivl Ice Shelf', 'Fimbul & Jelbart & Ekstrom Ice Shelves', 'Brunt & Riiser-Larsen Ice Shelves', 'Ross Ice Shelf'] @@ -69,7 +71,7 @@ def timeseries_massloss (file_path, log_path): # Read time data and convert from seconds to years id = Dataset(file_path, 'r') - new_time = id.variables['ocean_time'][:]/(365.25*24*60*60) + new_time = id.variables['ocean_time'][:]/(365.25*24*60*60) + add_years if exists(log_path): # Concatenate with time values from log file start_t = len(old_time) diff --git a/timeseries_seaice.py b/timeseries_seaice.py index 3196620..700eff7 100644 --- a/timeseries_seaice.py +++ b/timeseries_seaice.py @@ -11,7 +11,9 @@ # log_path = path to log file (if it exists, previously calculated values will # be read from it; regardless, it will be overwritten with all # calculated values following computation) -def timeseries_seaice (file_path, log_path): +# add_years = optional number of years to add to time array (multiple of 14 +# for repeating 1992-2005 spinup) +def timeseries_seaice (file_path, log_path, add_years=0): time = [] total_area = [] @@ -45,7 +47,7 @@ def timeseries_seaice (file_path, log_path): dx, dy = cartesian_grid_2d(lon, lat) dA = dx*dy # Read time values and convert from days to years - new_time = id.variables['time'][:]/365.25 + new_time = id.variables['time'][:]/365.25 + add_years # Concatenate with time values from log file for t in range(size(new_time)): time.append(new_time[t])