From b206bdcfafb5d502f690cc471e6a5a252aeeb5c6 Mon Sep 17 00:00:00 2001 From: Kaitlin Naughten Date: Fri, 8 Sep 2017 17:15:12 +1000 Subject: [PATCH] Updates to intercomparison figures --- common_grid.py | 8 +- depoorter_data | 33 +-- mip_aice_minmax_nsidc.py | 6 +- mip_dpt_calc.py | 8 +- mip_iceshelf_figures.py | 409 +++++++++++++++++++++++++++------ mip_mld.py | 123 +++++----- mip_regions_1var.py | 474 ++++++++++++++++++++------------------- mip_scatterplot.py | 3 +- 8 files changed, 696 insertions(+), 368 deletions(-) diff --git a/common_grid.py b/common_grid.py index d2a2978..cc23c91 100644 --- a/common_grid.py +++ b/common_grid.py @@ -284,12 +284,12 @@ def common_grid (roms_file, cice_file, out_file): # First calculate the two derivatives dsvstr_dx = ma.empty(shape(svstr_common)) # Forward difference approximation - dsvstr_dx[:,:-1] = (svstr_common[:,1:] - svstr_common[:,:-1])/(2*dx[:,:-1]) + dsvstr_dx[:,:-1] = (svstr_common[:,1:] - svstr_common[:,:-1])/dx[:,:-1] # Backward difference for the last row - dsvstr_dx[:,-1] = (svstr_common[:,-1] - svstr_common[:,-2])/(2*dx[:,-1]) + dsvstr_dx[:,-1] = (svstr_common[:,-1] - svstr_common[:,-2])/dx[:,-1] dsustr_dy = ma.empty(shape(sustr_common)) - dsustr_dy[:-1,:] = (sustr_common[1:,:] - sustr_common[:-1,:])/(2*dy[:-1,:]) - dsustr_dy[-1,:] = (sustr_common[-1,:] - sustr_common[-2,:])/(2*dy[-1,:]) + dsustr_dy[:-1,:] = (sustr_common[1:,:] - sustr_common[:-1,:])/dy[:-1,:] + dsustr_dy[-1,:] = (sustr_common[-1,:] - sustr_common[-2,:])/dy[-1,:] curl_str = dsvstr_dx - dsustr_dy curl_str = ma.masked_where(mask_common==0, curl_str) # Write to file diff --git a/depoorter_data b/depoorter_data index 5f0e1c4..c60acf8 100644 --- a/depoorter_data +++ b/depoorter_data @@ -1,21 +1,22 @@ -Name Mass loss Melt rate -Amery 39 +/- 21 0.65 +/- 0.35 -West 26 +/- 13 1.65 +/- 0.82 -Shackleton 76 +/- 23 2.20 +/- 0.67 +Name Mass loss Melt rate Rignot mass loss +Amery 39 +/- 21 0.65 +/- 0.35 35.5 +/- 23 +West 26 +/- 13 1.65 +/- 0.82 27.2 +/- 10 +Shackleton 76 +/- 23 2.20 +/- 0.67 72.6 +/- 15 Totten 64 +/- 12 9.89 +/- 1.92 Moscow Uni 28 +/- 7 4.93 +/- 1.32 -Mertz 5 +/- 4 0.87 +/- 0.79 -Ross 34 +/- 25 0.07 +/- 0.05 -Sulzberger 28 +/- 6 2.16 +/- 0.44 -Getz 136 +/- 23 4.09 +/- 0.68 -Thwaites 69 +/- 18 15.22 +/- 3.87 -Pine Island 95 +/- 14 15.96 +/- 2.38 -Abbot 86 +/- 22 2.72 +/- 0.7 -George VI 144 +/- 42 2.88 +/- 0.83 -Filchner-Ronne 50 +/- 40 0.12 +/- 0.09 -Brunt-RL 26 +/- 16 0.33 +/- 0.20 -Jelbart-Fimbul 24 +/- 13 0.52 +/- 0.27 -Total 1454 +/- 174 0.94 +/- 0.11 +Totten & MU 92 +/- 19 90.6 +/- 8 +Mertz 5 +/- 4 0.87 +/- 0.79 7.9 +/- 3 +Ross 34 +/- 25 0.07 +/- 0.05 47.7 +/- 34 +Sulzberger 28 +/- 6 2.16 +/- 0.44 18.2 +/- 3 +Getz 136 +/- 23 4.09 +/- 0.68 144.9 +/- 14 +Thwaites 69 +/- 18 15.22 +/- 3.87 97.5 +/- 7 +Pine Island 95 +/- 14 15.96 +/- 2.38 101.2 +/- 8 +Abbot 86 +/- 22 2.72 +/- 0.7 51.8 +/- 19 +George VI 144 +/- 42 2.88 +/- 0.83 135.4 +/- 40 (with Wilkins & Stange) +Filchner-Ronne 50 +/- 40 0.12 +/- 0.09 155.4 +/- 45 +Brunt-RL 26 +/- 16 0.33 +/- 0.20 9.7 +/- 16 +Jelbart-Fimbul 24 +/- 13 0.52 +/- 0.27 26.8 +/- 14 +Total 1454 +/- 174 0.94 +/- 0.11 Notes: Mertz a fair bit lower diff --git a/mip_aice_minmax_nsidc.py b/mip_aice_minmax_nsidc.py index 2e315e4..c36e8e2 100644 --- a/mip_aice_minmax_nsidc.py +++ b/mip_aice_minmax_nsidc.py @@ -360,14 +360,14 @@ def mip_aice_minmax_nsidc (cice_file, cice_log, fesom_mesh_path, fesom_output_di ax.set_yticks([]) # Add a colourbar at the bottom cbaxes = fig.add_axes([0.12, 0.04, 0.3, 0.04]) - cbar = colorbar(img, orientation='horizontal', ticks=arange(0,1+0.25,0.25), cax=cbaxes) + cbar = colorbar(img, orientation='horizontal', ticks=arange(0,1+0.25,0.25), cax=cbaxes, extend='min') cbar.ax.tick_params(labelsize=20) # Add extent timeseries on rightmost column, with more space for labels gs2 = GridSpec(2, 1) gs2.update(left=0.73, right=0.95, wspace=0.1, hspace=0.15) # February ax = subplot(gs2[0, 0]) - ax.plot(time_axis, nsidc_feb_extent, color='red', linewidth=1.5) + ax.plot(time_axis, nsidc_feb_extent, color='black', linewidth=2, linestyle='dashed') ax.plot(time_axis, cice_feb_extent, color='blue', linewidth=1.5) ax.plot(time_axis, fesom_feb_extent, color='green', linewidth=1.5) xlim([start_year, end_year]) @@ -379,7 +379,7 @@ def mip_aice_minmax_nsidc (cice_file, cice_log, fesom_mesh_path, fesom_output_di title('b) Sea ice extent\n'+r'(million km$^2$)', fontsize=26) # Extent timeseries, September ax = subplot(gs2[1, 0]) - ax.plot(time_axis, nsidc_sep_extent, color='red', label='NSIDC', linewidth=1.5) + ax.plot(time_axis, nsidc_sep_extent, color='black', label='NSIDC', linewidth=2, linestyle='dashed') ax.plot(time_axis, cice_sep_extent, color='blue', label='MetROMS', linewidth=1.5) ax.plot(time_axis, fesom_sep_extent, color='green', label='FESOM (high-res)', linewidth=1.5) xlim([start_year, end_year]) diff --git a/mip_dpt_calc.py b/mip_dpt_calc.py index 64e34d9..9d23c55 100644 --- a/mip_dpt_calc.py +++ b/mip_dpt_calc.py @@ -71,11 +71,11 @@ def mip_dpt_calc (roms_log, fesom_log_low, fesom_log_high): fesom_dpt_low = fesom_dpt_low[t_start_fesom:t_end_fesom] fesom_dpt_high = fesom_dpt_high[t_start_fesom:t_end_fesom] - # Calculate and print averages + # Calculate and print averages and standard deviations print 'Average Drake Passage Transport' - print 'MetROMS: ' + str(mean(roms_dpt)) - print 'FESOM low-res: ' + str(mean(fesom_dpt_low)) - print 'FESOM high-res: ' + str(mean(fesom_dpt_high)) + print 'MetROMS: ' + str(mean(roms_dpt)) + ' +/- ' + str(std(roms_dpt)) + print 'FESOM low-res: ' + str(mean(fesom_dpt_low)) + ' +/- ' + str(std(fesom_dpt_low)) + print 'FESOM high-res: ' + str(mean(fesom_dpt_high)) + ' +/- ' + str(std(fesom_dpt_high)) # Calculate and print trends # Also print p-values so we can see if it's statistically significant diff --git a/mip_iceshelf_figures.py b/mip_iceshelf_figures.py index e479a53..663fdbe 100644 --- a/mip_iceshelf_figures.py +++ b/mip_iceshelf_figures.py @@ -87,6 +87,155 @@ def get_min_max (roms_data, fesom_data_lr, fesom_data_hr, x_min, x_max, y_min, y return var_min, var_max +# Determine if the point (lon0, lat0) is within the triangular Element elm. +# Returns a boolean. +def in_triangle (elm, lon0, lat0): + + alpha = ((elm.lat[1] - elm.lat[2])*(lon0 - elm.lon[2]) + (elm.lon[2] - elm.lon[1])*(lat0 - elm.lat[2]))/((elm.lat[1] - elm.lat[2])*(elm.lon[0] - elm.lon[2]) + (elm.lon[2] - elm.lon[1])*(elm.lat[0] - elm.lat[2])) + beta = ((elm.lat[2] - elm.lat[0])*(lon0 - elm.lon[2]) + (elm.lon[0] - elm.lon[2])*(lat0 - elm.lat[2]))/((elm.lat[1] - elm.lat[2])*(elm.lon[0] - elm.lon[2]) + (elm.lon[2] - elm.lon[1])*(elm.lat[0] - elm.lat[2])) + gamma = 1 - alpha - beta + + return alpha >= 0 and beta >= 0 and gamma >= 0 + + +# Interpolate FESOM's bathymetry to a regular grid in the given region so that +# we can contour isobaths later. +# Input: +# lon_min, lon_max = bounds on longitude for regular grid (-180 to 180) +# lat_min, lat_max = bounds on latitude (-90 to 90) +# Output: +# x_reg, y_reg = polar coordinates of regular grid for plotting +# fesom_bathy_reg_lr, fesom_bathy_reg_hr = bathymetry interpolated to the +# regular grid, from the low-res and high-res mesh +# respectively +def fesom_interpolate_bathy (lon_min, lon_max, lat_min, lat_max): + + # Size of regular grid + num_lon = 500 + num_lat = 500 + + # Set up regular grid + # Start with boundaries + lon_reg_edges = linspace(lon_min, lon_max, num_lon+1) + lat_reg_edges = linspace(lat_min, lat_max, num_lat+1) + # Now get centres + lon_reg = 0.5*(lon_reg_edges[:-1] + lon_reg_edges[1:]) + lat_reg = 0.5*(lat_reg_edges[:-1] + lat_reg_edges[1:]) + # Make 2D versions + lon_reg_2d, lat_reg_2d = meshgrid(lon_reg, lat_reg) + # Calculate polar coordinates transformation for plotting + x_reg = -(lat_reg_2d+90)*cos(lon_reg_2d*deg2rad+pi/2) + y_reg = (lat_reg_2d+90)*sin(lon_reg_2d*deg2rad+pi/2) + # Set up array for bathymetry on this regular grid, for each FESOM mesh + fesom_bathy_reg_lr = zeros([num_lat, num_lon]) + fesom_bathy_reg_hr = zeros([num_lat, num_lon]) + + # Low-res + # For each element, check if a point on the regular lat-lon grid lies + # within. If so, do barycentric interpolation to that point, of bottom + # node depths (i.e. bathymetry). + for elm in elements_lr: + # Ignore ice shelf cavities + if not elm.cavity: + # Check if we are within domain of regular grid + if amax(elm.lon) > lon_min and amin(elm.lon) < lon_max and amax(elm.lat) > lat_min and amin(elm.lat) < lat_max: + # Find largest regular longitude value west of Element + tmp = nonzero(lon_reg > amin(elm.lon))[0] + if len(tmp) == 0: + # Element crosses the western boundary + iW = 0 + else: + iW = tmp[0] - 1 + # Find smallest regular longitude value east of Element + tmp = nonzero(lon_reg > amax(elm.lon))[0] + if len(tmp) == 0: + # Element crosses the eastern boundary + iE = num_lon + else: + iE = tmp[0] + # Find largest regular latitude value south of Element + tmp = nonzero(lat_reg > amin(elm.lat))[0] + if len(tmp) == 0: + # Element crosses the southern boundary + jS = 0 + else: + jS = tmp[0] - 1 + # Find smallest regular latitude value north of Element + tmp = nonzero(lat_reg > amax(elm.lat))[0] + if len(tmp) == 0: + # Element crosses the northern boundary + jN = num_lat + else: + jN = tmp[0] + for i in range(iW+1,iE): + for j in range(jS+1,jN): + # There is a chance that the regular gridpoint at (i,j) + # lies within this element + lon0 = lon_reg[i] + lat0 = lat_reg[j] + if in_triangle(elm, lon0, lat0): + # Yes it does + # Get area of entire triangle + area = triangle_area(elm.lon, elm.lat) + # Get area of each sub-triangle formed by + # (lon0, lat0) + area0 = triangle_area([lon0, elm.lon[1], elm.lon[2]], [lat0, elm.lat[1], elm.lat[2]]) + area1 = triangle_area([lon0, elm.lon[0], elm.lon[2]], [lat0, elm.lat[0], elm.lat[2]]) + area2 = triangle_area([lon0, elm.lon[0], elm.lon[1]], [lat0, elm.lat[0], elm.lat[1]]) + # Find fractional area of each + cff = [area0/area, area1/area, area2/area] + # Find bottom depth of each node + bathy_vals = [] + for n in range(3): + bathy_vals.append(elm.nodes[n].find_bottom().depth) + # Barycentric interpolation to lon0, lat0 + fesom_bathy_reg_lr[j,i] = sum(array(cff)*array(bathy_vals)) + # Mask out points which are identically zero (land and ice shelves) + fesom_bathy_reg_lr = ma.masked_where(fesom_bathy_reg_lr==0, fesom_bathy_reg_lr) + + # Repeat for high-res + for elm in elements_hr: + if not elm.cavity: + if amax(elm.lon) > lon_min and amin(elm.lon) < lon_max and amax(elm.lat) > lat_min and amin(elm.lat) < lat_max: + tmp = nonzero(lon_reg > amin(elm.lon))[0] + if len(tmp) == 0: + iW = 0 + else: + iW = tmp[0] - 1 + tmp = nonzero(lon_reg > amax(elm.lon))[0] + if len(tmp) == 0: + iE = num_lon + else: + iE = tmp[0] + tmp = nonzero(lat_reg > amin(elm.lat))[0] + if len(tmp) == 0: + jS = 0 + else: + jS = tmp[0] - 1 + tmp = nonzero(lat_reg > amax(elm.lat))[0] + if len(tmp) == 0: + jN = num_lat + else: + jN = tmp[0] + for i in range(iW+1,iE): + for j in range(jS+1,jN): + lon0 = lon_reg[i] + lat0 = lat_reg[j] + if in_triangle(elm, lon0, lat0): + area = triangle_area(elm.lon, elm.lat) + area0 = triangle_area([lon0, elm.lon[1], elm.lon[2]], [lat0, elm.lat[1], elm.lat[2]]) + area1 = triangle_area([lon0, elm.lon[0], elm.lon[2]], [lat0, elm.lat[0], elm.lat[2]]) + area2 = triangle_area([lon0, elm.lon[0], elm.lon[1]], [lat0, elm.lat[0], elm.lat[1]]) + cff = [area0/area, area1/area, area2/area] + bathy_vals = [] + for n in range(3): + bathy_vals.append(elm.nodes[n].find_bottom().depth) + fesom_bathy_reg_hr[j,i] = sum(array(cff)*array(bathy_vals)) + fesom_bathy_reg_hr = ma.masked_where(fesom_bathy_reg_hr==0, fesom_bathy_reg_hr) + + return x_reg, y_reg, fesom_bathy_reg_lr, fesom_bathy_reg_hr + + # Create a 2D vector field for vertically averaged velocity in MetROMS, low-res # FESOM, and high-res FESOM for the given region. Average velocity over # horizontal bins (in x and y) for easy plotting. @@ -259,7 +408,79 @@ def plot_draft (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter): ylim([y_min, y_max]) ax.set_xticks([]) ax.set_yticks([]) - # Colourbar on the right + # Colourbar + cbar = colorbar(img, cax=cbaxes, ticks=cbar_ticks) + + +# Plot water column thickness for MetROMS, low-res FESOM, and high-res FESOM +# for the given region, in the given axes. +# Input: +# x_min, x_max, y_min, y_max = bounds on x and y (using polar coordinate +# transformation from above) for the desired region +# gs = GridSpec object of size 1x3 to plot in +# cbaxes = Axes object for location of colourbar +# cbar_ticks = 1D array containing values for ticks on colourbar +# letter = 'a', 'b', 'c', etc. to add before the variable title, for +# use in a figure showing multiple variables +def plot_wct (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter): + + # Set up a grey square for FESOM to fill the background with land + x_reg_fesom, y_reg_fesom = meshgrid(linspace(x_min, x_max, num=100), linspace(y_min, y_max, num=100)) + land_square = zeros(shape(x_reg_fesom)) + # Find bounds on variable in this region + var_min, var_max = get_min_max(roms_wct, fesom_wct_lr, fesom_wct_hr, x_min, x_max, y_min, y_max) + + # MetROMS + ax = subplot(gs[0,0], aspect='equal') + # First shade land and zice in grey + grey_cmap = ListedColormap([(0.6, 0.6, 0.6)]) + pcolor(roms_x, roms_y, land_zice, cmap=grey_cmap) + # Fill in the missing circle + pcolor(x_reg_roms, y_reg_roms, land_circle, cmap=grey_cmap) + # Now shade the data + pcolor(roms_x, roms_y, roms_wct, vmin=var_min, vmax=var_max, cmap='jet') + xlim([x_min, x_max]) + ylim([y_min, y_max]) + ax.set_xticks([]) + ax.set_yticks([]) + + # FESOM low-res + ax = subplot(gs[0,1], aspect='equal') + # Start with land background + contourf(x_reg_fesom, y_reg_fesom, land_square, 1, colors=(('0.6', '0.6', '0.6'))) + # Add ice shelf elements + img = PatchCollection(patches_lr, cmap='jet') + img.set_array(array(fesom_wct_lr)) + img.set_edgecolor('face') + img.set_clim(vmin=var_min, vmax=var_max) + ax.add_collection(img) + # Mask out the open ocean in white + overlay = PatchCollection(mask_patches_lr, facecolor=(1,1,1)) + overlay.set_edgecolor('face') + ax.add_collection(overlay) + xlim([x_min, x_max]) + ylim([y_min, y_max]) + ax.set_xticks([]) + ax.set_yticks([]) + # Main title + title(letter + ') Water column thickness (m)', fontsize=20) + + # FESOM high-res + ax = subplot(gs[0,2], aspect='equal') + contourf(x_reg_fesom, y_reg_fesom, land_square, 1, colors=(('0.6', '0.6', '0.6'))) + img = PatchCollection(patches_hr, cmap='jet') + img.set_array(array(fesom_wct_hr)) + img.set_edgecolor('face') + img.set_clim(vmin=var_min, vmax=var_max) + ax.add_collection(img) + overlay = PatchCollection(mask_patches_hr, facecolor=(1,1,1)) + overlay.set_edgecolor('face') + ax.add_collection(overlay) + xlim([x_min, x_max]) + ylim([y_min, y_max]) + ax.set_xticks([]) + ax.set_yticks([]) + # Colourbar cbar = colorbar(img, cax=cbaxes, ticks=cbar_ticks) @@ -375,7 +596,7 @@ def plot_melt (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, change_points ylim([y_min, y_max]) ax.set_xticks([]) ax.set_yticks([]) - # Colourbar on the right + # Colourbar cbar = colorbar(img, cax=cbaxes, ticks=cbar_ticks) text(0.5, y0, 'FESOM (high-res)', fontsize=18, horizontalalignment='center', transform=ax.transAxes) @@ -390,7 +611,9 @@ def plot_melt (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, change_points # cbar_ticks = 1D array containing values for ticks on colourbar # letter = 'a', 'b', 'c', etc. to add before the bottom water temp title, for # use in a figure showing multiple variables -def plot_bwtemp (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter): +# bathy_contour = optional float containing single isobath to contour (positive, +# in metres; make sure you call fesom_interpolate_bathy first) +def plot_bwtemp (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter, bathy_contour=None): # Set up a grey square for FESOM to fill the background with land x_reg_fesom, y_reg_fesom = meshgrid(linspace(x_min, x_max, num=100), linspace(y_min, y_max, num=100)) @@ -410,6 +633,9 @@ def plot_bwtemp (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter): # Add a black contour for the ice shelf front rcParams['contour.negative_linestyle'] = 'solid' contour(roms_x, roms_y, zice_contour, levels=[min_zice], colors=('black')) + if bathy_contour is not None: + # Dashed line to contour isobath + contour(roms_x, roms_y, roms_bathy, levels=[bathy_contour], colors=('black'), linestyles=('dashed')) xlim([x_min, x_max]) ylim([y_min, y_max]) ax.set_xticks([]) @@ -428,6 +654,9 @@ def plot_bwtemp (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter): # Add ice shelf front contour lines fesom_contours_lr = LineCollection(contour_lines_lr, edgecolor='black', linewidth=1) ax.add_collection(fesom_contours_lr) + if bathy_contour is not None: + # Overlay contour on regular grid + contour(x_reg, y_reg, fesom_bathy_reg_lr, levels=[bathy_contour], colors=('black'), linestyles=('dashed')) xlim([x_min, x_max]) ylim([y_min, y_max]) ax.set_xticks([]) @@ -443,16 +672,15 @@ def plot_bwtemp (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter): img.set_edgecolor('face') img.set_clim(vmin=var_min, vmax=var_max) ax.add_collection(img) - #overlay = PatchCollection(mask_patches_hr, facecolor=(1,1,1)) - #overlay.set_edgecolor('face') - #ax.add_collection(overlay) fesom_contours_hr = LineCollection(contour_lines_hr, edgecolor='black', linewidth=1) ax.add_collection(fesom_contours_hr) + if bathy_contour is not None: + contour(x_reg, y_reg, fesom_bathy_reg_hr, levels=[bathy_contour], colors=('black'), linestyles=('dashed')) xlim([x_min, x_max]) ylim([y_min, y_max]) ax.set_xticks([]) ax.set_yticks([]) - # Colourbar on the right + # Colourbar cbar = colorbar(img, cax=cbaxes, ticks=cbar_ticks) @@ -466,7 +694,9 @@ def plot_bwtemp (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter): # cbar_ticks = 1D array containing values for ticks on colourbar # letter = 'a', 'b', 'c', etc. to add before the bottom water salinity title, # for use in a figure showing multiple variables -def plot_bwsalt (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter): +# bathy_contour = optional float containing single isobath to contour (positive, +# in metres; make sure you call fesom_interpolate_bathy first) +def plot_bwsalt (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter, bathy_contour=None): # Set up a grey square for FESOM to fill the background with land x_reg_fesom, y_reg_fesom = meshgrid(linspace(x_min, x_max, num=100), linspace(y_min, y_max, num=100)) @@ -486,6 +716,9 @@ def plot_bwsalt (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter): # Add a black contour for the ice shelf front rcParams['contour.negative_linestyle'] = 'solid' contour(roms_x, roms_y, zice_contour, levels=[min_zice], colors=('black')) + if bathy_contour is not None: + # Dashed line to contour isobath + contour(roms_x, roms_y, roms_bathy, levels=[bathy_contour], colors=('black'), linestyles=('dashed')) xlim([x_min, x_max]) ylim([y_min, y_max]) ax.set_xticks([]) @@ -501,13 +734,12 @@ def plot_bwsalt (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter): img.set_edgecolor('face') img.set_clim(vmin=var_min, vmax=var_max) ax.add_collection(img) - # Mask out the open ocean in white - #overlay = PatchCollection(mask_patches_lr, facecolor=(1,1,1)) - #overlay.set_edgecolor('face') - #ax.add_collection(overlay) # Add ice shelf front contour lines fesom_contours_lr = LineCollection(contour_lines_lr, edgecolor='black', linewidth=1) ax.add_collection(fesom_contours_lr) + if bathy_contour is not None: + # Overlay contour on regular grid + contour(x_reg, y_reg, fesom_bathy_reg_lr, levels=[bathy_contour], colors=('black'), linestyles=('dashed')) xlim([x_min, x_max]) ylim([y_min, y_max]) ax.set_xticks([]) @@ -523,16 +755,15 @@ def plot_bwsalt (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter): img.set_edgecolor('face') img.set_clim(vmin=var_min, vmax=var_max) ax.add_collection(img) - #overlay = PatchCollection(mask_patches_hr, facecolor=(1,1,1)) - #overlay.set_edgecolor('face') - #ax.add_collection(overlay) fesom_contours_hr = LineCollection(contour_lines_hr, edgecolor='black', linewidth=1) ax.add_collection(fesom_contours_hr) xlim([x_min, x_max]) ylim([y_min, y_max]) + if bathy_contour is not None: + contour(x_reg, y_reg, fesom_bathy_reg_hr, levels=[bathy_contour], colors=('black'), linestyles=('dashed')) ax.set_xticks([]) ax.set_yticks([]) - # Colourbar on the right + # Colourbar cbar = colorbar(img, cax=cbaxes, ticks=cbar_ticks) @@ -554,7 +785,11 @@ def plot_bwsalt (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter): # variables (eg "George VI Ice Shelf" for the Bellingshausen plot) # arrow_scale, arrow_headwidth, arrow_headlength = optional parameters for # arrows on vector overlay -def plot_velavg (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr, letter, loc_string=None, arrow_scale=0.9, arrow_headwidth=8, arrow_headlength=9): +# model_titles = optional boolean indicating to add model titles to the top +# of the plot (only for Filchner-Ronne). Default False. +# y0 = if model_titles=True, y-coordinate of model titles. Play around between +# 1.15 and 1.35. Default 1.25. +def plot_velavg (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr, letter, loc_string=None, arrow_scale=0.9, arrow_headwidth=8, arrow_headlength=9, model_titles=False, y0=1.25): # Set up a grey square for FESOM to fill the background with land x_reg_fesom, y_reg_fesom = meshgrid(linspace(x_min, x_max, num=100), linspace(y_min, y_max, num=100)) @@ -571,12 +806,15 @@ def plot_velavg (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, x_centres, pcolor(x_reg_roms, y_reg_roms, land_circle, cmap=grey_cmap) # Now shade the data pcolor(roms_x, roms_y, roms_speed, vmin=var_min, vmax=var_max, cmap='cool') - # Overlay vectors - quiver(x_centres, y_centres, roms_ubin, roms_vbin, scale=arrow_scale, headwidth=arrow_headwidth, headlength=arrow_headlength, color='black') + if x_centres is not None: + # Overlay vectors + 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]) ax.set_xticks([]) ax.set_yticks([]) + if model_titles: + text(0.5, y0, 'MetROMS', fontsize=18, horizontalalignment='center', transform=ax.transAxes) # FESOM low-res ax = subplot(gs[0,1], aspect='equal') @@ -592,12 +830,15 @@ def plot_velavg (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, x_centres, overlay = PatchCollection(mask_patches_lr, facecolor=(1,1,1)) overlay.set_edgecolor('face') ax.add_collection(overlay) - # Overlay vectors - quiver(x_centres, y_centres, fesom_ubin_lr, fesom_vbin_lr, scale=arrow_scale, headwidth=arrow_headwidth, headlength=arrow_headlength, color='black') + if x_centres is not None: + # Overlay vectors + 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]) ax.set_xticks([]) ax.set_yticks([]) + if model_titles: + text(0.5, y0, 'FESOM (low-res)', fontsize=18, horizontalalignment='center', transform=ax.transAxes) # Main title if loc_string is None: title(letter + ') Vertically averaged ocean velocity (m/s)', fontsize=20) @@ -615,13 +856,16 @@ def plot_velavg (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, x_centres, overlay = PatchCollection(mask_patches_hr, facecolor=(1,1,1)) overlay.set_edgecolor('face') ax.add_collection(overlay) - quiver(x_centres, y_centres, fesom_ubin_hr, fesom_vbin_hr, scale=arrow_scale, headwidth=arrow_headwidth, headlength=arrow_headlength, color='black') + if x_centres is not None: + 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]) ax.set_xticks([]) - ax.set_yticks([]) - # Colourbar on the right + ax.set_yticks([]) + # Colourbar cbar = colorbar(img, cax=cbaxes, ticks=cbar_ticks) + if model_titles: + text(0.5, y0, 'FESOM (high-res)', fontsize=18, horizontalalignment='center', transform=ax.transAxes) # Plot zonal slices (latitude vs depth) of temperature and salinity for @@ -849,6 +1093,8 @@ def plot_zonal_ts (lon0, lat_min, lat_max, lat_ticks, lat_labels, depth_min, dep open_ocn = copy(roms_mask) open_ocn[roms_zice!=0] = 0 land_zice = ma.masked_where(open_ocn==1, open_ocn) +# Mask land and zice out of bathymetry array for isobath contours +roms_bathy = ma.masked_where(open_ocn==0, roms_h) # Convert grid to spherical coordinates roms_x = -(roms_lat+90)*cos(roms_lon*deg2rad+pi/2) roms_y = (roms_lat+90)*sin(roms_lon*deg2rad+pi/2) @@ -928,7 +1174,7 @@ def plot_zonal_ts (lon0, lat_min, lat_max, lat_ticks, lat_labels, depth_min, dep # FESOM low-res # Calculate draft at each element, averaged over 3 corners -# Equivalent to depth of surfce layer +# Equivalent to depth of surface layer fesom_draft_lr = [] for elm in elements_lr: if elm.cavity: @@ -940,6 +1186,28 @@ def plot_zonal_ts (lon0, lat_min, lat_max, lat_ticks, lat_labels, depth_min, dep if elm.cavity: fesom_draft_hr.append(mean([elm.nodes[0].depth, elm.nodes[1].depth, elm.nodes[2].depth])) +print 'Calculating water column thickness' + +# ROMS +# Add h (positive) and zice (negative) +roms_wct = roms_h + roms_zice +# Mask the open ocean and land +roms_wct = ma.masked_where(roms_zice==0, roms_wct) + +# FESOM low-res +# Calculate wct at each element, averaged over 3 corners +# Equivalent to depth of bottom layer minus depth of surface layer +fesom_wct_lr = [] +for elm in elements_lr: + if elm.cavity: + fesom_wct_lr.append(mean([elm.nodes[0].find_bottom().depth - elm.nodes[0].depth, elm.nodes[1].find_bottom().depth - elm.nodes[1].depth, elm.nodes[2].find_bottom().depth - elm.nodes[2].depth])) + +# FESOM high-res +fesom_wct_hr = [] +for elm in elements_hr: + if elm.cavity: + fesom_wct_hr.append(mean([elm.nodes[0].find_bottom().depth - elm.nodes[0].depth, elm.nodes[1].find_bottom().depth - elm.nodes[1].depth, elm.nodes[2].find_bottom().depth - elm.nodes[2].depth])) + print 'Calculating ice shelf melt rate' # ROMS @@ -1249,39 +1517,45 @@ def plot_zonal_ts (lon0, lat_min, lat_max, lat_ticks, lat_labels, depth_min, dep x_max_tmp = -4.5 y_min_tmp = 1 y_max_tmp = 10 -fig = figure(figsize=(8,14)) +fig = figure(figsize=(8,22)) fig.patch.set_facecolor('white') -# Melt rate +# Melt gs_a = GridSpec(1,3) -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]) +gs_a.update(left=0.05, right=0.9, bottom=0.764, top=0.89, wspace=0.05) +cbaxes_tmp = fig.add_axes([0.91, 0.777, 0.025, 0.1]) 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, [-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) -gs_b.update(left=0.05, right=0.9, bottom=0.555, top=0.71, wspace=0.05) -cbaxes_tmp = fig.add_axes([0.91, 0.575, 0.025, 0.12]) +gs_b.update(left=0.05, right=0.9, bottom=0.614, top=0.74, wspace=0.05) +cbaxes_tmp = fig.add_axes([0.91, 0.627, 0.025, 0.11]) cbar_ticks = arange(0, 0.2+0.1, 0.1) 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.9, 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.375, top=0.53, wspace=0.05) -cbaxes_tmp = fig.add_axes([0.91, 0.395, 0.025, 0.12]) +gs_c.update(left=0.05, right=0.9, bottom=0.464, top=0.59, wspace=0.05) +cbaxes_tmp = fig.add_axes([0.91, 0.477, 0.025, 0.1]) 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) -gs_d.update(left=0.05, right=0.9, bottom=0.195, top=0.35, wspace=0.05) -cbaxes_tmp = fig.add_axes([0.91, 0.215, 0.025, 0.12]) +gs_d.update(left=0.05, right=0.9, bottom=0.314, top=0.44, wspace=0.05) +cbaxes_tmp = fig.add_axes([0.91, 0.327, 0.025, 0.1]) cbar_ticks = arange(34.3, 34.7+0.2, 0.2) plot_bwsalt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_d, cbaxes_tmp, cbar_ticks, 'd') -# Ice shelf draft +# Water column thickness gs_e = GridSpec(1,3) -gs_e.update(left=0.05, right=0.9, bottom=0.015, top=0.17, wspace=0.05) -cbaxes_tmp = fig.add_axes([0.91, 0.035, 0.025, 0.12]) +gs_e.update(left=0.05, right=0.9, bottom=0.164, top=0.29, wspace=0.05) +cbaxes_tmp = fig.add_axes([0.91, 0.177, 0.025, 0.1]) +cbar_ticks = arange(200, 1000+400, 400) +plot_wct(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_e, cbaxes_tmp, cbar_ticks, 'e') +# Ice shelf draft +gs_f = GridSpec(1,3) +gs_f.update(left=0.05, right=0.9, bottom=0.014, top=0.14, wspace=0.05) +cbaxes_tmp = fig.add_axes([0.91, 0.027, 0.025, 0.1]) cbar_ticks = arange(500, 1500+500, 500) -plot_draft(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_e, cbaxes_tmp, cbar_ticks, 'e') +plot_draft(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_f, cbaxes_tmp, cbar_ticks, 'f') suptitle('Filchner-Ronne Ice Shelf', fontsize=30) fig.show() fig.savefig('filchner_ronne.png') @@ -1300,21 +1574,21 @@ def plot_zonal_ts (lon0, lat_min, lat_max, lat_ticks, lat_labels, depth_min, dep 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, [-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) +#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) gs_b.update(left=0.11, right=0.9, bottom=0.57, top=0.72, wspace=0.05) cbaxes_tmp = fig.add_axes([0.91, 0.6, 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.9, 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, None, None, None, None, None, None, None, None, 'b', arrow_scale=0.9, arrow_headwidth=8, arrow_headlength=9) # Temp and salt slices through 1W: Fimbul Ice Shelf lat_min = -71.5 lat_max = -69.5 lat_ticks = [-71.5, -71, -70.5, -70, -69.5] lat_labels = ['', r'71$^{\circ}$S', '', r'70$^{\circ}$S', ''] -depth_min = -2400 +depth_min = -1500 depth_max = 0 -depth_ticks = [-2000, -1500, -1000, -500, 0] -depth_labels = ['2000', '1500', '1000', '500', '0'] +depth_ticks = [-1500, -1000, -500, 0] +depth_labels = ['1500', '1000', '500', '0'] gs_c = GridSpec(1,3) gs_c.update(left=0.11, right=0.9, bottom=0.33, top=0.53, wspace=0.05) cbaxes_tmp1 = fig.add_axes([0.91, 0.36, 0.025, 0.14]) @@ -1333,39 +1607,45 @@ def plot_zonal_ts (lon0, lat_min, lat_max, lat_ticks, lat_labels, depth_min, dep x_max_tmp = 20.5 y_min_tmp = 4.75 y_max_tmp = 8 -fig = figure(figsize=(10,12)) +fig = figure(figsize=(10,14)) fig.patch.set_facecolor('white') # Melt gs_a = GridSpec(1,3) -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]) +gs_a.update(left=0.05, right=0.9, bottom=0.764, top=0.89, wspace=0.05) +cbaxes_tmp = fig.add_axes([0.91, 0.777, 0.025, 0.1]) cbar_ticks = arange(0, 40+20, 20) plot_melt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_a, cbaxes_tmp, cbar_ticks, [2, 5, 11], 'a', 1.25, [70], [-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, 20, 15) gs_b = GridSpec(1,3) -gs_b.update(left=0.05, right=0.9, bottom=0.555, top=0.71, wspace=0.05) -cbaxes_tmp = fig.add_axes([0.91, 0.575, 0.025, 0.12]) +gs_b.update(left=0.05, right=0.9, bottom=0.614, top=0.74, wspace=0.05) +cbaxes_tmp = fig.add_axes([0.91, 0.627, 0.025, 0.11]) cbar_ticks = arange(0.1, 0.3+0.1, 0.1) 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.9, arrow_headwidth=7, arrow_headlength=8) # Bottom water temperature 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]) +gs_c.update(left=0.05, right=0.9, bottom=0.464, top=0.59, wspace=0.05) +cbaxes_tmp = fig.add_axes([0.91, 0.477, 0.025, 0.1]) cbar_ticks = arange(-2.4, -1.4+0.5, 0.5) 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) -gs_d.update(left=0.05, right=0.9, bottom=0.195, top=0.35, wspace=0.05) -cbaxes_tmp = fig.add_axes([0.91, 0.215, 0.025, 0.12]) +gs_d.update(left=0.05, right=0.9, bottom=0.314, top=0.44, wspace=0.05) +cbaxes_tmp = fig.add_axes([0.91, 0.327, 0.025, 0.1]) cbar_ticks = arange(34, 34.4+0.2, 0.2) plot_bwsalt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_d, cbaxes_tmp, cbar_ticks, 'd') -# Ice shelf draft +# Water column thickness gs_e = GridSpec(1,3) -gs_e.update(left=0.05, right=0.9, bottom=0.015, top=0.17, wspace=0.05) -cbaxes_tmp = fig.add_axes([0.91, 0.035, 0.025, 0.12]) +gs_e.update(left=0.05, right=0.9, bottom=0.164, top=0.29, wspace=0.05) +cbaxes_tmp = fig.add_axes([0.91, 0.177, 0.025, 0.1]) +cbar_ticks = arange(250, 1250+500, 500) +plot_wct(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_e, cbaxes_tmp, cbar_ticks, 'e') +# Ice shelf draft +gs_f = GridSpec(1,3) +gs_f.update(left=0.05, right=0.9, bottom=0.014, top=0.14, wspace=0.05) +cbaxes_tmp = fig.add_axes([0.91, 0.027, 0.025, 0.1]) cbar_ticks = arange(200, 2200+1000, 1000) -plot_draft(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_e, cbaxes_tmp, cbar_ticks, 'e') +plot_draft(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_f, cbaxes_tmp, cbar_ticks, 'f') suptitle('Amery Ice Shelf', fontsize=30) fig.show() fig.savefig('amery.png') @@ -1384,17 +1664,18 @@ def plot_zonal_ts (lon0, lat_min, lat_max, lat_ticks, lat_labels, depth_min, dep 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.15, [90, 120], [-65]) # Bottom water temperature +x_reg, y_reg, fesom_bathy_reg_lr, fesom_bathy_reg_hr = fesom_interpolate_bathy(75, 160, -70, -63) 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.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') +plot_bwtemp(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_b, cbaxes_tmp, cbar_ticks, 'b', bathy_contour=1000) # 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.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') +plot_bwsalt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_c, cbaxes_tmp, cbar_ticks, 'c', bathy_contour=1000) # Velocity for Totten x_min_tmp = 20 x_max_tmp = 21.5 @@ -1472,17 +1753,18 @@ def plot_zonal_ts (lon0, lat_min, lat_max, lat_ticks, lat_labels, depth_min, dep cbar_ticks = arange(0, 14+7, 7) plot_melt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_a, cbaxes_tmp, cbar_ticks, [1, 3, 6], 'a', 1.15, [-130, -110], [-75]) # Bottom water temperature +x_reg, y_reg, fesom_bathy_reg_lr, fesom_bathy_reg_hr = fesom_interpolate_bathy(-140, -90, -76, -72) 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.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') +plot_bwtemp(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_b, cbaxes_tmp, cbar_ticks, 'b', bathy_contour=1000) # 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.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') +plot_bwsalt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_c, cbaxes_tmp, cbar_ticks, 'c', bathy_contour=1000) # Velocity for PIG x_min_tmp = -15.6 x_max_tmp = -14.1 @@ -1512,17 +1794,18 @@ def plot_zonal_ts (lon0, lat_min, lat_max, lat_ticks, lat_labels, depth_min, dep cbar_ticks = arange(0, 12+6, 6) 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.15, [-100, -80], [-70]) # Bottom water temperature +x_reg, y_reg, fesom_bathy_reg_lr, fesom_bathy_reg_hr = fesom_interpolate_bathy(-110, -60, -76, -66) 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.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') +plot_bwtemp(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_b, cbaxes_tmp, cbar_ticks, 'b', bathy_contour=1000) # 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.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') +plot_bwsalt(x_min_tmp, x_max_tmp, y_min_tmp, y_max_tmp, gs_c, cbaxes_tmp, cbar_ticks, 'c', bathy_contour=1000) # Velocity for George VI x_min_tmp = -18.75 x_max_tmp = -15.5 diff --git a/mip_mld.py b/mip_mld.py index 4219ef0..3d36295 100644 --- a/mip_mld.py +++ b/mip_mld.py @@ -150,58 +150,82 @@ def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file fesom_mld[season,:] = array(mld_season) print 'Plotting' - fig = figure(figsize=(19,9)) - # MetROMS, summer, ACC - ax = fig.add_subplot(2, 4, 1, aspect='equal') + # ACC + fig1 = figure(figsize=(9.5,9)) + # Summer + # MetROMS + ax = fig1.add_subplot(2, 2, 1, aspect='equal') pcolor(roms_x, roms_y, roms_mld[0,:,:], vmin=0, vmax=max_bound_summer, cmap='jet') - text(-62, 0, 'MetROMS', fontsize=24, ha='right') - title(season_names[0], fontsize=24) + text(-65, 0, season_names[0], fontsize=24, ha='right') + title('MetROMS', fontsize=24) xlim([-nbdry1, nbdry1]) ylim([-nbdry1, nbdry1]) # 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) + text(x_ticks[i], y_ticks[i], lon_labels[i], ha='center', rotation=lon_rot[i], fontsize=12) 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]) - 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) + # FESOM + ax = fig1.add_subplot(2, 2, 2, aspect='equal') + img = PatchCollection(patches, cmap='jet') + img.set_array(fesom_mld[0,:]) + img.set_clim(vmin=0, vmax=max_bound_summer) + img.set_edgecolor('face') + ax.add_collection(img) xlim([-nbdry1, nbdry1]) ylim([-nbdry1, nbdry1]) ax.set_xticks([]) ax.set_yticks([]) - # MetROMS, winter, continental shelf - ax = fig.add_subplot(2, 4, 4, aspect='equal') + title('FESOM (high-res)', fontsize=24) + # Add a colorbar for summer + cbaxes = fig1.add_axes([0.91, 0.55, 0.02, 0.3]) + cbar = colorbar(img, cax=cbaxes, extend='max', ticks=arange(0, max_bound_summer+50, 50)) + cbar.ax.tick_params(labelsize=20) + # Winter + # MetROMS + ax = fig1.add_subplot(2, 2, 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([-nbdry2, nbdry2]) - ylim([-nbdry2, nbdry2]) + text(-65, 0, season_names[2], fontsize=24, ha='right') + xlim([-nbdry1, nbdry1]) + ylim([-nbdry1, nbdry1]) ax.set_xticks([]) ax.set_yticks([]) - # FESOM, summer, ACC - ax = fig.add_subplot(2, 4, 5, aspect='equal') + # FESOM + ax = fig1.add_subplot(2, 2, 4, aspect='equal') img = PatchCollection(patches, cmap='jet') - img.set_array(fesom_mld[0,:]) - img.set_clim(vmin=0, vmax=max_bound_summer) + img.set_array(fesom_mld[2,:]) + img.set_clim(vmin=0, vmax=max_bound_winter) img.set_edgecolor('face') ax.add_collection(img) xlim([-nbdry1, nbdry1]) ylim([-nbdry1, nbdry1]) 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 - ax = fig.add_subplot(2, 4, 6, aspect='equal') + # Add a colorbar for winter + cbaxes = fig1.add_axes([0.91, 0.15, 0.02, 0.3]) + cbar = colorbar(img, cax=cbaxes, extend='max', ticks=arange(0, max_bound_winter+200, 200)) + cbar.ax.tick_params(labelsize=20) + # Add the main title + suptitle('Mixed layer depth (m), 2002-2016 average', fontsize=30) + # Decrease space between plots + subplots_adjust(wspace=0.025,hspace=0.025) + fig1.show() + fig1.savefig('mld_acc.png') + + # Continental shelf + fig2 = figure(figsize=(9.5,9)) + # Summer + # MetROMS + ax = fig2.add_subplot(2, 2, 1, aspect='equal') + pcolor(roms_x, roms_y, roms_mld[0,:,:], vmin=0, vmax=max_bound_summer, cmap='jet') + text(-28, 0, season_names[0], fontsize=24, ha='right') + title('MetROMS', fontsize=24) + xlim([-nbdry2, nbdry2]) + ylim([-nbdry2, nbdry2]) + ax.set_xticks([]) + ax.set_yticks([]) + # FESOM + ax = fig2.add_subplot(2, 2, 2, aspect='equal') img = PatchCollection(patches, cmap='jet') img.set_array(fesom_mld[0,:]) img.set_clim(vmin=0, vmax=max_bound_summer) @@ -211,23 +235,22 @@ def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file ylim([-nbdry2, nbdry2]) 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)) + title('FESOM (high-res)', fontsize=24) + # Add a colorbar for summer + cbaxes = fig2.add_axes([0.91, 0.55, 0.02, 0.3]) + cbar = colorbar(img, cax=cbaxes, extend='max', ticks=arange(0, max_bound_summer+50, 50)) cbar.ax.tick_params(labelsize=20) - # FESOM, winter, ACC - ax = fig.add_subplot(2, 4, 7, aspect='equal') - img = PatchCollection(patches, cmap='jet') - img.set_array(fesom_mld[2,:]) - img.set_clim(vmin=0, vmax=max_bound_winter) - img.set_edgecolor('face') - ax.add_collection(img) - xlim([-nbdry1, nbdry1]) - ylim([-nbdry1, nbdry1]) + # Winter + # MetROMS + ax = fig2.add_subplot(2, 2, 3, aspect='equal') + pcolor(roms_x, roms_y, roms_mld[2,:,:], vmin=0, vmax=max_bound_winter, cmap='jet') + text(-28, 0, season_names[2], fontsize=24, ha='right') + xlim([-nbdry2, nbdry2]) + ylim([-nbdry2, nbdry2]) ax.set_xticks([]) ax.set_yticks([]) - # FESOM, winter, continental shelf - ax = fig.add_subplot(2, 4, 8, aspect='equal') + # FESOM + ax = fig2.add_subplot(2, 2, 4, aspect='equal') img = PatchCollection(patches, cmap='jet') img.set_array(fesom_mld[2,:]) img.set_clim(vmin=0, vmax=max_bound_winter) @@ -237,16 +260,16 @@ def mip_mld (roms_grid, roms_seasonal_file, fesom_mesh_path, fesom_seasonal_file ylim([-nbdry2, nbdry2]) 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)) + # Add a colorbar for winter + cbaxes = fig2.add_axes([0.91, 0.15, 0.02, 0.3]) + cbar = colorbar(img, cax=cbaxes, extend='max', ticks=arange(0, max_bound_winter+200, 200)) cbar.ax.tick_params(labelsize=20) # Add the main title suptitle('Mixed layer depth (m), 2002-2016 average', fontsize=30) # Decrease space between plots subplots_adjust(wspace=0.025,hspace=0.025) - #fig.show() - fig.savefig('mld.png') + fig2.show() + fig2.savefig('mld_shelf.png') # Command-line interface diff --git a/mip_regions_1var.py b/mip_regions_1var.py index 9f2387b..4bef5f0 100644 --- a/mip_regions_1var.py +++ b/mip_regions_1var.py @@ -16,12 +16,13 @@ from unrotate_grid import * # For each of the 8 ice shelf regions analysed in this intercomparison paper, -# make one 3x1 figure for 5 different variables (ice shelf draft, ice shelf +# make one 3x1 figure for 7 different variables (ice shelf draft, ice shelf # melt rate, bottom water temperature, bottom water salinity, vertically -# averaged velocity) where each variable is averaged over the years 2002-2016 -# and zoomed into the region of interest. The 3 parts of each figure are -# MetROMS, FESOM low-res, and FESOM high-res. The velocity figures show -# absolute value (speed) in colours, and direction with vector arrows. +# averaged velocity, bathymetry, water column thickness) where each variable is +# averaged over the years 2002-2016 and zoomed into the region of interest. The +# 3 parts of each figure are MetROMS, FESOM low-res, and FESOM high-res. The +# velocity figures show absolute value (speed) in colours, and direction with +# vector arrows. def mip_regions_1var (): # Path to ROMS grid file @@ -51,7 +52,7 @@ def mip_regions_1var (): # Size of each plot in the y direction ysize = [8, 6, 7, 9, 7, 9, 10, 7] # Variables to process - var_names = ['vel'] #['draft', 'melt', 'temp', 'salt', 'vel'] + var_names = ['draft', 'melt', 'temp', 'salt', 'vel', 'bathy', 'wct'] # Constants sec_per_year = 365*24*3600 deg2rad = pi/180.0 @@ -110,6 +111,12 @@ def mip_regions_1var (): if var == 'draft': # Swap sign on existing zice field; nothing more to read roms_data = -1*roms_zice + elif var == 'bathy': + # Point to h field; nothing more to read + roms_data = roms_h + elif var == 'wct': + # Add h (positive) and zice (negative); nothing more to read + roms_data = roms_h + roms_zice else: id = Dataset(roms_file, 'r') if var == 'melt': @@ -164,11 +171,9 @@ def mip_regions_1var (): # Get speed roms_data = sqrt(u_rho**2 + v_rho**2) id.close() - # Mask the open ocean and land out of the data field - roms_data = ma.masked_where(roms_zice==0, roms_data) print 'Reading FESOM low-res fields' - if var != 'draft': + if var not in ['draft', 'bathy', 'wct']: if var == 'melt': id = Dataset(fesom_file_lr_i, 'r') # Convert from m/s to m/y @@ -275,6 +280,13 @@ def mip_regions_1var (): if var == 'draft': # Ice shelf draft is depth of surface layer fesom_data_lr.append(mean([elm.nodes[0].depth, elm.nodes[1].depth, elm.nodes[2].depth])) + elif var == 'bathy': + # Bathymetry is depth of bottom layer + fesom_data_lr.append(mean([elm.nodes[0].find_bottom().depth, elm.nodes[1].find_bottom().depth, elm.nodes[2].find_bottom().depth])) + elif var == 'wct': + # Water column thickness is depth of bottom layer minus + # depth of surface layer + fesom_data_lr.append(mean([elm.nodes[0].find_bottom().depth - elm.nodes[0].depth, elm.nodes[1].find_bottom().depth - elm.nodes[1].depth, elm.nodes[2].find_bottom().depth - elm.nodes[2].depth])) elif var in ['melt', 'vel']: # Surface nodes fesom_data_lr.append(mean([node_data_lr[elm.nodes[0].id], node_data_lr[elm.nodes[1].id], node_data_lr[elm.nodes[2].id]])) @@ -284,7 +296,7 @@ def mip_regions_1var (): print 'Reading FESOM high-res fields' # As before - if var != 'draft': + if var not in ['draft', 'bathy', 'wct']: if var == 'melt': id = Dataset(fesom_file_hr_i, 'r') node_data_hr = id.variables['wnet'][0,:]*sec_per_year @@ -365,231 +377,239 @@ def mip_regions_1var (): if elm.cavity: if var == 'draft': fesom_data_hr.append(mean([elm.nodes[0].depth, elm.nodes[1].depth, elm.nodes[2].depth])) + elif var == 'bathy': + fesom_data_hr.append(mean([elm.nodes[0].find_bottom().depth, elm.nodes[1].find_bottom().depth, elm.nodes[2].find_bottom().depth])) + elif var == 'wct': + fesom_data_hr.append(mean([elm.nodes[0].find_bottom().depth - elm.nodes[0].depth, elm.nodes[1].find_bottom().depth - elm.nodes[1].depth, elm.nodes[2].find_bottom().depth - elm.nodes[2].depth])) elif var in ['melt', 'vel']: fesom_data_hr.append(mean([node_data_hr[elm.nodes[0].id], node_data_hr[elm.nodes[1].id], node_data_hr[elm.nodes[2].id]])) elif var in ['temp', 'salt']: fesom_data_hr.append(mean([node_data_hr[elm.nodes[0].find_bottom().id], node_data_hr[elm.nodes[1].find_bottom().id], node_data_hr[elm.nodes[2].find_bottom().id]])) - # Loop over regions - for index in range(num_regions): - print 'Processing ' + region_names[index] - # Set up a grey square for FESOM to fill the background with land - x_reg_fesom, y_reg_fesom = meshgrid(linspace(x_min[index], x_max[index], num=100), linspace(y_min[index], y_max[index], num=100)) - land_square = zeros(shape(x_reg_fesom)) - # Find bounds on variable in this region, for both ROMS and FESOM - # Start with ROMS - loc = (roms_x >= x_min[index])*(roms_x <= x_max[index])*(roms_y >= y_min[index])*(roms_y <= y_max[index]) - var_min = amin(roms_data[loc]) - var_max = amax(roms_data[loc]) - # Modify with FESOM - # Low-res - i = 0 - for elm in elements_lr: - if elm.cavity: - if any(elm.x >= x_min[index]) and any(elm.x <= x_max[index]) and any(elm.y >= y_min[index]) and any(elm.y <= y_max[index]): - if fesom_data_lr[i] < var_min: - var_min = fesom_data_lr[i] - if fesom_data_lr[i] > var_max: - var_max = fesom_data_lr[i] - i += 1 - # High-res - i = 0 - for elm in elements_hr: - if elm.cavity: - if any(elm.x >= x_min[index]) and any(elm.x <= x_max[index]) and any(elm.y >= y_min[index]) and any(elm.y <= y_max[index]): - if fesom_data_hr[i] < var_min: - var_min = fesom_data_hr[i] - if fesom_data_hr[i] > var_max: - var_max = fesom_data_hr[i] - i += 1 - if var == 'melt': - # Special colour map - if var_min < 0: - # There is refreezing here; include blue for elements < 0 - cmap_vals = array([var_min, 0, 0.25*var_max, 0.5*var_max, 0.75*var_max, var_max]) - cmap_colors = [(0.26, 0.45, 0.86), (1, 1, 1), (1, 0.9, 0.4), (0.99, 0.59, 0.18), (0.5, 0.0, 0.08), (0.96, 0.17, 0.89)] - cmap_vals_norm = (cmap_vals - var_min)/(var_max - var_min) - cmap_list = [] - for i in range(size(cmap_vals)): - cmap_list.append((cmap_vals_norm[i], cmap_colors[i])) - mf_cmap = LinearSegmentedColormap.from_list('melt_freeze', cmap_list) + # Loop over regions + for index in range(num_regions): + print 'Processing ' + region_names[index] + # Set up a grey square for FESOM to fill the background with land + x_reg_fesom, y_reg_fesom = meshgrid(linspace(x_min[index], x_max[index], num=100), linspace(y_min[index], y_max[index], num=100)) + land_square = zeros(shape(x_reg_fesom)) + # Find bounds on variable in this region, for both ROMS and FESOM + # Start with ROMS + loc = (roms_x >= x_min[index])*(roms_x <= x_max[index])*(roms_y >= y_min[index])*(roms_y <= y_max[index]) + var_min = amin(roms_data[loc]) + var_max = amax(roms_data[loc]) + # Modify with FESOM + # Low-res + i = 0 + for elm in elements_lr: + if elm.cavity: + if any(elm.x >= x_min[index]) and any(elm.x <= x_max[index]) and any(elm.y >= y_min[index]) and any(elm.y <= y_max[index]): + if fesom_data_lr[i] < var_min: + var_min = fesom_data_lr[i] + if fesom_data_lr[i] > var_max: + var_max = fesom_data_lr[i] + i += 1 + # High-res + i = 0 + for elm in elements_hr: + if elm.cavity: + if any(elm.x >= x_min[index]) and any(elm.x <= x_max[index]) and any(elm.y >= y_min[index]) and any(elm.y <= y_max[index]): + if fesom_data_hr[i] < var_min: + var_min = fesom_data_hr[i] + if fesom_data_hr[i] > var_max: + var_max = fesom_data_hr[i] + i += 1 + if var == 'melt': + # Special colour map + if var_min < 0: + # There is refreezing here; include blue for elements < 0 + cmap_vals = array([var_min, 0, 0.25*var_max, 0.5*var_max, 0.75*var_max, var_max]) + cmap_colors = [(0.26, 0.45, 0.86), (1, 1, 1), (1, 0.9, 0.4), (0.99, 0.59, 0.18), (0.5, 0.0, 0.08), (0.96, 0.17, 0.89)] + cmap_vals_norm = (cmap_vals - var_min)/(var_max - var_min) + cmap_list = [] + for i in range(size(cmap_vals)): + cmap_list.append((cmap_vals_norm[i], cmap_colors[i])) + mf_cmap = LinearSegmentedColormap.from_list('melt_freeze', cmap_list) + else: + # No refreezing + cmap_vals = array([0, 0.25*var_max, 0.5*var_max, 0.75*var_max, var_max]) + cmap_colors = [(1, 1, 1), (1, 0.9, 0.4), (0.99, 0.59, 0.18), (0.5, 0.0, 0.08), (0.96, 0.17, 0.89)] + cmap_vals_norm = cmap_vals/var_max + cmap_list = [] + for i in range(size(cmap_vals)): + cmap_list.append((cmap_vals_norm[i], cmap_colors[i])) + mf_cmap = LinearSegmentedColormap.from_list('melt_freeze', cmap_list) + colour_map = mf_cmap + elif var == 'vel': + colour_map = 'cool' else: - # No refreezing - cmap_vals = array([0, 0.25*var_max, 0.5*var_max, 0.75*var_max, var_max]) - cmap_colors = [(1, 1, 1), (1, 0.9, 0.4), (0.99, 0.59, 0.18), (0.5, 0.0, 0.08), (0.96, 0.17, 0.89)] - cmap_vals_norm = cmap_vals/var_max - cmap_list = [] - for i in range(size(cmap_vals)): - cmap_list.append((cmap_vals_norm[i], cmap_colors[i])) - mf_cmap = LinearSegmentedColormap.from_list('melt_freeze', cmap_list) - colour_map = mf_cmap - elif var == 'vel': - colour_map = 'cool' - else: - colour_map = 'jet' - if var == 'vel': - # Make vectors for overlay - # Set up bins (edges) - x_bins = linspace(x_min[index], x_max[index], num=num_bins+1) - y_bins = linspace(y_min[index], y_max[index], num=num_bins+1) - # Calculate centres of bins (for plotting) - x_centres = 0.5*(x_bins[:-1] + x_bins[1:]) - y_centres = 0.5*(y_bins[:-1] + y_bins[1:]) - # ROMS - # First set up arrays to integrate velocity in each bin - # Simple averaging of all the points inside each bin - roms_u = zeros([size(y_centres), size(x_centres)]) - roms_v = zeros([size(y_centres), size(x_centres)]) - roms_num_pts = zeros([size(y_centres), size(x_centres)]) - # First convert to polar coordinates, rotate to account for - # longitude in circumpolar projection, and convert back to vector - # components - theta_roms = arctan2(v_rho, u_rho) - theta_circ_roms = theta_roms - roms_lon*deg2rad - u_circ_roms = roms_data*cos(theta_circ_roms) # roms_data is speed - v_circ_roms = roms_data*sin(theta_circ_roms) - # Loop over all points (can't find a better way to do this) - for j in range(size(roms_data,0)): - for i in range(size(roms_data,1)): - # Make sure data isn't masked (i.e. land or open ocean) - if u_circ_roms[j,i] is not ma.masked: - # Check if we're in the region of interest - if roms_x[j,i] > x_min[index] and roms_x[j,i] < x_max[index] and roms_y[j,i] > y_min[index] and roms_y[j,i] < y_max[index]: - # Figure out which bins this falls into - x_index = nonzero(x_bins > roms_x[j,i])[0][0]-1 - y_index = nonzero(y_bins > roms_y[j,i])[0][0]-1 - # Integrate - roms_u[y_index, x_index] += u_circ_roms[j,i] - roms_v[y_index, x_index] += v_circ_roms[j,i] - roms_num_pts[y_index, x_index] += 1 - # Convert from sums to averages - # First mask out points with no data - roms_u = ma.masked_where(roms_num_pts==0, roms_u) - roms_v = ma.masked_where(roms_num_pts==0, roms_v) - # Divide everything else by the number of points - flag = roms_num_pts > 0 - roms_u[flag] = roms_u[flag]/roms_num_pts[flag] - roms_v[flag] = roms_v[flag]/roms_num_pts[flag] + colour_map = 'jet' + if var == 'vel': + # Make vectors for overlay + # Set up bins (edges) + x_bins = linspace(x_min[index], x_max[index], num=num_bins+1) + y_bins = linspace(y_min[index], y_max[index], num=num_bins+1) + # Calculate centres of bins (for plotting) + x_centres = 0.5*(x_bins[:-1] + x_bins[1:]) + y_centres = 0.5*(y_bins[:-1] + y_bins[1:]) + # ROMS + # First set up arrays to integrate velocity in each bin + # Simple averaging of all the points inside each bin + roms_u = zeros([size(y_centres), size(x_centres)]) + roms_v = zeros([size(y_centres), size(x_centres)]) + roms_num_pts = zeros([size(y_centres), size(x_centres)]) + # First convert to polar coordinates, rotate to account for + # longitude in circumpolar projection, and convert back to vector + # components + theta_roms = arctan2(v_rho, u_rho) + theta_circ_roms = theta_roms - roms_lon*deg2rad + u_circ_roms = roms_data*cos(theta_circ_roms) # roms_data is speed + v_circ_roms = roms_data*sin(theta_circ_roms) + # Loop over all points (can't find a better way to do this) + for j in range(size(roms_data,0)): + for i in range(size(roms_data,1)): + # Make sure data isn't masked (i.e. land or open ocean) + if u_circ_roms[j,i] is not ma.masked: + # Check if we're in the region of interest + if roms_x[j,i] > x_min[index] and roms_x[j,i] < x_max[index] and roms_y[j,i] > y_min[index] and roms_y[j,i] < y_max[index]: + # Figure out which bins this falls into + x_index = nonzero(x_bins > roms_x[j,i])[0][0]-1 + y_index = nonzero(y_bins > roms_y[j,i])[0][0]-1 + # Integrate + roms_u[y_index, x_index] += u_circ_roms[j,i] + roms_v[y_index, x_index] += v_circ_roms[j,i] + roms_num_pts[y_index, x_index] += 1 + # Convert from sums to averages + # First mask out points with no data + roms_u = ma.masked_where(roms_num_pts==0, roms_u) + roms_v = ma.masked_where(roms_num_pts==0, roms_v) + # Divide everything else by the number of points + flag = roms_num_pts > 0 + roms_u[flag] = roms_u[flag]/roms_num_pts[flag] + roms_v[flag] = roms_v[flag]/roms_num_pts[flag] + # FESOM low-res + fesom_u_lr = zeros([size(y_centres), size(x_centres)]) + fesom_v_lr = zeros([size(y_centres), size(x_centres)]) + fesom_num_pts_lr = zeros([size(y_centres), size(x_centres)]) + theta_fesom_lr = arctan2(node_v_lr, node_u_lr) + theta_circ_fesom_lr = theta_fesom_lr - fesom_lon_lr*deg2rad + u_circ_fesom_lr = node_data_lr*cos(theta_circ_fesom_lr) # node_data is speed + v_circ_fesom_lr = node_data_lr*sin(theta_circ_fesom_lr) + # Loop over 2D nodes to fill in the velocity bins as before + for n in range(fesom_n2d_lr): + if fesom_cavity_lr[n]: + if fesom_x_lr[n] > x_min[index] and fesom_x_lr[n] < x_max[index] and fesom_y_lr[n] > y_min[index] and fesom_y_lr[n] < y_max[index]: + x_index = nonzero(x_bins > fesom_x_lr[n])[0][0]-1 + y_index = nonzero(y_bins > fesom_y_lr[n])[0][0]-1 + fesom_u_lr[y_index, x_index] += u_circ_fesom_lr[n] + fesom_v_lr[y_index, x_index] += v_circ_fesom_lr[n] + fesom_num_pts_lr[y_index, x_index] += 1 + fesom_u_lr = ma.masked_where(fesom_num_pts_lr==0, fesom_u_lr) + fesom_v_lr = ma.masked_where(fesom_num_pts_lr==0, fesom_v_lr) + flag = fesom_num_pts_lr > 0 + fesom_u_lr[flag] = fesom_u_lr[flag]/fesom_num_pts_lr[flag] + fesom_v_lr[flag] = fesom_v_lr[flag]/fesom_num_pts_lr[flag] + # FESOM high-res + fesom_u_hr = zeros([size(y_centres), size(x_centres)]) + fesom_v_hr = zeros([size(y_centres), size(x_centres)]) + fesom_num_pts_hr = zeros([size(y_centres), size(x_centres)]) + theta_fesom_hr = arctan2(node_v_hr, node_u_hr) + theta_circ_fesom_hr = theta_fesom_hr - fesom_lon_hr*deg2rad + u_circ_fesom_hr = node_data_hr*cos(theta_circ_fesom_hr) # node_data is speed + v_circ_fesom_hr = node_data_hr*sin(theta_circ_fesom_hr) + # Loop over 2D nodes to fill in the velocity bins as before + for n in range(fesom_n2d_hr): + if fesom_cavity_hr[n]: + if fesom_x_hr[n] > x_min[index] and fesom_x_hr[n] < x_max[index] and fesom_y_hr[n] > y_min[index] and fesom_y_hr[n] < y_max[index]: + x_index = nonzero(x_bins > fesom_x_hr[n])[0][0]-1 + y_index = nonzero(y_bins > fesom_y_hr[n])[0][0]-1 + fesom_u_hr[y_index, x_index] += u_circ_fesom_hr[n] + fesom_v_hr[y_index, x_index] += v_circ_fesom_hr[n] + fesom_num_pts_hr[y_index, x_index] += 1 + fesom_u_hr = ma.masked_where(fesom_num_pts_hr==0, fesom_u_hr) + fesom_v_hr = ma.masked_where(fesom_num_pts_hr==0, fesom_v_hr) + flag = fesom_num_pts_hr > 0 + fesom_u_hr[flag] = fesom_u_hr[flag]/fesom_num_pts_hr[flag] + fesom_v_hr[flag] = fesom_v_hr[flag]/fesom_num_pts_hr[flag] + # Plot + fig = figure(figsize=(20, ysize[index])) + fig.patch.set_facecolor('white') + # MetROMS + ax = fig.add_subplot(1,3,1, aspect='equal') + # First shade land and zice in grey + contourf(roms_x, roms_y, land_zice, 1, colors=(('0.6', '0.6', '0.6'))) + # Fill in the missing circle + contourf(x_reg_roms, y_reg_roms, land_circle, 1, colors=(('0.6', '0.6', '0.6'))) + # Now shade the data + pcolor(roms_x, roms_y, roms_data, vmin=var_min, vmax=var_max, cmap=colour_map) + if var == 'vel': + # Overlay vectors + quiver(x_centres, y_centres, roms_u, roms_v, scale=1.5, headwidth=6, headlength=7, color='black') + xlim([x_min[index], x_max[index]]) + ylim([y_min[index], y_max[index]]) + axis('off') + title('MetROMS', fontsize=24) # FESOM low-res - fesom_u_lr = zeros([size(y_centres), size(x_centres)]) - fesom_v_lr = zeros([size(y_centres), size(x_centres)]) - fesom_num_pts_lr = zeros([size(y_centres), size(x_centres)]) - theta_fesom_lr = arctan2(node_v_lr, node_u_lr) - theta_circ_fesom_lr = theta_fesom_lr - fesom_lon_lr*deg2rad - u_circ_fesom_lr = node_data_lr*cos(theta_circ_fesom_lr) # node_data is speed - v_circ_fesom_lr = node_data_lr*sin(theta_circ_fesom_lr) - # Loop over 2D nodes to fill in the velocity bins as before - for n in range(fesom_n2d_lr): - if fesom_cavity_lr[n]: - if fesom_x_lr[n] > x_min[index] and fesom_x_lr[n] < x_max[index] and fesom_y_lr[n] > y_min[index] and fesom_y_lr[n] < y_max[index]: - x_index = nonzero(x_bins > fesom_x_lr[n])[0][0]-1 - y_index = nonzero(y_bins > fesom_y_lr[n])[0][0]-1 - fesom_u_lr[y_index, x_index] += u_circ_fesom_lr[n] - fesom_v_lr[y_index, x_index] += v_circ_fesom_lr[n] - fesom_num_pts_lr[y_index, x_index] += 1 - fesom_u_lr = ma.masked_where(fesom_num_pts_lr==0, fesom_u_lr) - fesom_v_lr = ma.masked_where(fesom_num_pts_lr==0, fesom_v_lr) - flag = fesom_num_pts_lr > 0 - fesom_u_lr[flag] = fesom_u_lr[flag]/fesom_num_pts_lr[flag] - fesom_v_lr[flag] = fesom_v_lr[flag]/fesom_num_pts_lr[flag] + ax = fig.add_subplot(1,3,2, aspect='equal') + # Start with land background + contourf(x_reg_fesom, y_reg_fesom, land_square, 1, colors=(('0.6', '0.6', '0.6'))) + # Add ice shelf elements + img = PatchCollection(patches_lr, cmap=colour_map) + img.set_array(array(fesom_data_lr)) + img.set_edgecolor('face') + img.set_clim(vmin=var_min, vmax=var_max) + ax.add_collection(img) + # Mask out the open ocean in white + overlay = PatchCollection(mask_patches_lr, facecolor=(1,1,1)) + overlay.set_edgecolor('face') + ax.add_collection(overlay) + if var == 'vel': + # Overlay vectors + quiver(x_centres, y_centres, fesom_u_lr, fesom_v_lr, scale=1.5, headwidth=6, headlength=7, color='black') + xlim([x_min[index], x_max[index]]) + ylim([y_min[index], y_max[index]]) + axis('off') + title('FESOM (low-res)', fontsize=24) # FESOM high-res - fesom_u_hr = zeros([size(y_centres), size(x_centres)]) - fesom_v_hr = zeros([size(y_centres), size(x_centres)]) - fesom_num_pts_hr = zeros([size(y_centres), size(x_centres)]) - theta_fesom_hr = arctan2(node_v_hr, node_u_hr) - theta_circ_fesom_hr = theta_fesom_hr - fesom_lon_hr*deg2rad - u_circ_fesom_hr = node_data_hr*cos(theta_circ_fesom_hr) # node_data is speed - v_circ_fesom_hr = node_data_hr*sin(theta_circ_fesom_hr) - # Loop over 2D nodes to fill in the velocity bins as before - for n in range(fesom_n2d_hr): - if fesom_cavity_hr[n]: - if fesom_x_hr[n] > x_min[index] and fesom_x_hr[n] < x_max[index] and fesom_y_hr[n] > y_min[index] and fesom_y_hr[n] < y_max[index]: - x_index = nonzero(x_bins > fesom_x_hr[n])[0][0]-1 - y_index = nonzero(y_bins > fesom_y_hr[n])[0][0]-1 - fesom_u_hr[y_index, x_index] += u_circ_fesom_hr[n] - fesom_v_hr[y_index, x_index] += v_circ_fesom_hr[n] - fesom_num_pts_hr[y_index, x_index] += 1 - fesom_u_hr = ma.masked_where(fesom_num_pts_hr==0, fesom_u_hr) - fesom_v_hr = ma.masked_where(fesom_num_pts_hr==0, fesom_v_hr) - flag = fesom_num_pts_hr > 0 - fesom_u_hr[flag] = fesom_u_hr[flag]/fesom_num_pts_hr[flag] - fesom_v_hr[flag] = fesom_v_hr[flag]/fesom_num_pts_hr[flag] - # Plot - fig = figure(figsize=(20, ysize[index])) - fig.patch.set_facecolor('white') - # MetROMS - ax = fig.add_subplot(1,3,1, aspect='equal') - # First shade land and zice in grey - contourf(roms_x, roms_y, land_zice, 1, colors=(('0.6', '0.6', '0.6'))) - # Fill in the missing circle - contourf(x_reg_roms, y_reg_roms, land_circle, 1, colors=(('0.6', '0.6', '0.6'))) - # Now shade the data - pcolor(roms_x, roms_y, roms_data, vmin=var_min, vmax=var_max, cmap=colour_map) - if var == 'vel': - # Overlay vectors - quiver(x_centres, y_centres, roms_u, roms_v, scale=1.5, headwidth=6, headlength=7, color='black') - xlim([x_min[index], x_max[index]]) - ylim([y_min[index], y_max[index]]) - axis('off') - title('MetROMS', fontsize=24) - # FESOM low-res - ax = fig.add_subplot(1,3,2, aspect='equal') - # Start with land background - contourf(x_reg_fesom, y_reg_fesom, land_square, 1, colors=(('0.6', '0.6', '0.6'))) - # Add ice shelf elements - img = PatchCollection(patches_lr, cmap=colour_map) - img.set_array(array(fesom_data_lr)) - img.set_edgecolor('face') - img.set_clim(vmin=var_min, vmax=var_max) - ax.add_collection(img) - # Mask out the open ocean in white - overlay = PatchCollection(mask_patches_lr, facecolor=(1,1,1)) - overlay.set_edgecolor('face') - ax.add_collection(overlay) - if var == 'vel': - # Overlay vectors - quiver(x_centres, y_centres, fesom_u_lr, fesom_v_lr, scale=1.5, headwidth=6, headlength=7, color='black') - xlim([x_min[index], x_max[index]]) - ylim([y_min[index], y_max[index]]) - axis('off') - title('FESOM (low-res)', fontsize=24) - # FESOM high-res - ax = fig.add_subplot(1,3,3, aspect='equal') - contourf(x_reg_fesom, y_reg_fesom, land_square, 1, colors=(('0.6', '0.6', '0.6'))) - img = PatchCollection(patches_hr, cmap=colour_map) - img.set_array(array(fesom_data_hr)) - img.set_edgecolor('face') - img.set_clim(vmin=var_min, vmax=var_max) - ax.add_collection(img) - overlay = PatchCollection(mask_patches_hr, facecolor=(1,1,1)) - overlay.set_edgecolor('face') - ax.add_collection(overlay) - if var == 'vel': - # Overlay vectors - quiver(x_centres, y_centres, fesom_u_hr, fesom_v_hr, scale=1.5, headwidth=6, headlength=7, color='black') - xlim([x_min[index], x_max[index]]) - ylim([y_min[index], y_max[index]]) - axis('off') - title('FESOM (high-res)', fontsize=24) - # Colourbar on the right - cbaxes = fig.add_axes([0.92, 0.2, 0.01, 0.6]) - cbar = colorbar(img, cax=cbaxes) - cbar.ax.tick_params(labelsize=20) - # Main title - if var == 'draft': - title_string = ' draft (m)' - elif var == 'melt': - title_string = ' melt rate (m/y)' - elif var == 'temp': - title_string = r' bottom water temperature ($^{\circ}$C)' - elif var == 'salt': - title_string = ' bottom water salinity (psu)' - elif var == 'vel': - title_string = ' vertically averaged ocean velocity (m/s)' - suptitle(region_names[index] + title_string, fontsize=30) - subplots_adjust(wspace=0.05) - #fig.show() - fig.savefig(fig_heads[index] + '_' + var + '.png') + ax = fig.add_subplot(1,3,3, aspect='equal') + contourf(x_reg_fesom, y_reg_fesom, land_square, 1, colors=(('0.6', '0.6', '0.6'))) + img = PatchCollection(patches_hr, cmap=colour_map) + img.set_array(array(fesom_data_hr)) + img.set_edgecolor('face') + img.set_clim(vmin=var_min, vmax=var_max) + ax.add_collection(img) + overlay = PatchCollection(mask_patches_hr, facecolor=(1,1,1)) + overlay.set_edgecolor('face') + ax.add_collection(overlay) + if var == 'vel': + # Overlay vectors + quiver(x_centres, y_centres, fesom_u_hr, fesom_v_hr, scale=1.5, headwidth=6, headlength=7, color='black') + xlim([x_min[index], x_max[index]]) + ylim([y_min[index], y_max[index]]) + axis('off') + title('FESOM (high-res)', fontsize=24) + # Colourbar on the right + cbaxes = fig.add_axes([0.92, 0.2, 0.01, 0.6]) + cbar = colorbar(img, cax=cbaxes) + cbar.ax.tick_params(labelsize=20) + # Main title + if var == 'draft': + title_string = ' draft (m)' + elif var == 'bathy': + title_string = ' bathymetry (m)' + elif var == 'wct': + title_string = ' water column thickness (m)' + elif var == 'melt': + title_string = ' melt rate (m/y)' + elif var == 'temp': + title_string = r' bottom water temperature ($^{\circ}$C)' + elif var == 'salt': + title_string = ' bottom water salinity (psu)' + elif var == 'vel': + title_string = ' vertically averaged ocean velocity (m/s)' + suptitle(region_names[index] + title_string, fontsize=30) + subplots_adjust(wspace=0.05) + #fig.show() + fig.savefig(fig_heads[index] + '_' + var + '.png') # Command-line interface diff --git a/mip_scatterplot.py b/mip_scatterplot.py index 0a456ff..3ce6f32 100644 --- a/mip_scatterplot.py +++ b/mip_scatterplot.py @@ -11,7 +11,7 @@ def mip_scatterplot (roms_logfile, fesom_logfile_lr, fesom_logfile_hr): # Number of output steps per year in FESOM peryear = 365/5 # Name of each ice shelf - names = ['Larsen D', 'Larsen C', 'Wilkins & George VI & Stange', 'Filchner-Ronne', 'Abbot', 'Pine Island Glacier', 'Thwaites', 'Dotson', 'Getz', 'Nickerson', 'Sulzberger', 'Mertz', 'Totten & Moscow University', 'Shackleton', 'West', 'Amery', 'Prince Harald', 'Baudouin & Borchgrevink', 'Lazarev', 'Nivl', 'Fimbul & Jelbart & Ekstrom', 'Brunt & Riiser-Larsen', 'Ross'] + names = ['Larsen D', 'Larsen C', 'Wilkins & George VI & Stange', 'Filchner-Ronne', 'Abbot', 'Pine Island', 'Thwaites', 'Dotson', 'Getz', 'Nickerson', 'Sulzberger', 'Mertz', 'Totten & Moscow University', 'Shackleton', 'West', 'Amery', 'Prince Harald', 'Baudouin & Borchgrevink', 'Lazarev', 'Nivl', 'Fimbul & Jelbart & Ekstrom', 'Brunt & Riiser-Larsen', 'Ross'] # Observed mass loss (Rignot 2013) and uncertainty for each ice shelf, in Gt/y obs_massloss = [1.4, 20.7, 135.4, 155.4, 51.8, 101.2, 97.5, 45.2, 144.9, 4.2, 18.2, 7.9, 90.6, 72.6, 27.2, 35.5, -2, 21.6, 6.3, 3.9, 26.8, 9.7, 47.7] obs_massloss_error = [14, 67, 40, 45, 19, 8, 7, 4, 14, 2, 3, 3, 8, 15, 10, 23, 3, 18, 2, 2, 14, 16, 34] @@ -177,6 +177,7 @@ def mip_scatterplot (roms_logfile, fesom_logfile_lr, fesom_logfile_hr): legend(numpoints=1,loc='lower left') setp(gca().get_legend().get_texts(), fontsize='13') fig.show() + fig.savefig('scatterplot.png') # Command-line interface