diff --git a/file_guide.txt b/file_guide.txt index 0a52653..f64ea7f 100644 --- a/file_guide.txt +++ b/file_guide.txt @@ -1379,12 +1379,12 @@ mip_calc_massloss.py: Calculate the average mass loss for each ice shelf over mip_regions_1var.py: 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 melt rate, bottom - water temperature, bottom water salinity, surface 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 surface velocity figures show absolute value (speed) - in colours, and direction with vector arrows. + 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. To run: First clone my "fesomtools" GitHub repository and replace '/short/y99/kaa561/fesomtools' (near the top of the file) with the path to the cloned @@ -1394,6 +1394,48 @@ mip_regions_1var.py: For each of the 8 ice shelf regions analysed in this and type "run mip_regions_1var.py". The script will output a bunch of pngs. +mip_iceshelf_figures.py: This is the giant monster script to generate 8 + multi-part figures showing ice shelf processes in the + 8 regions defined in the intercomparison paper. Each + figure is size Nx3 where N=3,4,5 is the number of + variables shown (current options are ice shelf melt + rate, ice shelf draft, bottom water temperature or + salinity, vertically averaged velocity, zonal slices + of temperature and salinity) and 3 is the number of + models (MetROMS, low-res FESOM, high-res FESOM). All + variables are averaged over the period 2002-2016, i.e. + the first 10 years of the simulation are discarded as + spinup, and seasonal cycles and interannual + variability are not considered. + In order to avoid unnecessarily repetitive code, this + script contains multiple functions which accept + location bounds as arguments (eg to calculate min/max + over the given region, to plot ice shelf draft over + the given region for each model, etc.) but in order to + prevent the function argument lists from getting out + of control, most of the variables are global and so + the majority of this script is not wrapped within a + function. Be very careful if you are running this + script in an active python session to make sure there + are no conflicting variable names previously defined. + To run: First clone my "fesomtools" GitHub repository + and replace '/short/y99/kaa561/fesomtools' + (near the top of the file) with the path to + the cloned repository on your system. Also + make sure the paths to the ROMS and FESOM + files near the top of the script are correct. + To use the predetermined regions and variable + combinations, simply open python or ipython + and type "run mip_iceshelf_figures.py". The + script will output a bunch of pngs. To modify + for other regions or variable combinations, + scroll down to "USER MODIFIED SECTION" to see + examples of how to call the functions to + create multi-part figures. You will have to + fiddle around with the x and y bounds, + GridSpec arguments, colourbar axes, etc. so + that the positioning looks okay. + ***UTILITY FUNCTIONS*** diff --git a/mip_iceshelf_figures.py b/mip_iceshelf_figures.py index 772119f..d57422c 100644 --- a/mip_iceshelf_figures.py +++ b/mip_iceshelf_figures.py @@ -19,6 +19,31 @@ from unrotate_vector import * from unrotate_grid import * +# This is the giant monster script to generate 8 multi-part figures showing +# ice shelf processes in the 8 regions defined in the intercomparison paper. +# Each figure is size Nx3 where N=3,4,5 is the number of variables shown +# (current options are ice shelf melt rate, ice shelf draft, bottom water +# temperature or salinity, vertically averaged velocity, zonal slices of +# temperature and salinity) and 3 is the number of models (MetROMS, low-res +# FESOM, high-res FESOM). All variables are averaged over the period 2002-2016, +# i.e. the first 10 years of the simulation are discarded as spinup, and +# seasonal cycles and interannual variability are not considered. + +# In order to avoid unnecessarily repetitive code, this script contains +# multiple functions which accept location bounds as arguments (eg to calculate +# min/max over the given region, to plot ice shelf draft over the given region +# for each model, etc.) but in order to prevent the function argument lists from +# getting out of control, most of the variables are global and so the majority +# of this script is not wrapped within a function. Be very careful if you are +# running this script in an active python session to make sure there are no +# conflicting variable names previously defined. + +# To modify for your purposes, adjust the file paths (immediately below) as +# needed, and then scroll down to "USER MODIFIED SECTION" to see examples of +# how to call the functions to create multi-part figures. You will have to +# fiddle around with the x and y bounds, GridSpec arguments, colourbar axes, +# etc. so that the positioning looks okay. + # File paths roms_grid = '/short/m68/kaa561/metroms_iceshelf/apps/common/grid/circ30S_quarterdegree.nc' roms_file = '/short/m68/kaa561/metroms_iceshelf/tmproms/run/intercomparison/2002_2016_avg.nc' @@ -727,6 +752,16 @@ fig.savefig('larsen.png') +# For the given variable in MetROMS, low-res FESOM, and high-res FESOM, find +# the max and min values across all 3 models in the given region. +# Input: +# roms_data = 2D field of data on the MetROMS grid (lat x lon) +# fesom_data_lr, fesom_data_hr = data for each 2D element on the FESOM low-res +# and high-res meshes respectively +# x_min, x_max, y_min, y_max = bounds on x and y (using polar coordinate +# transformation from above) for the desired region +# Output: +# var_min, var_max = min and max data values in this region across all 3 models def get_min_max (roms_data, fesom_data_lr, fesom_data_hr, x_min, x_max, y_min, y_max): # Start with ROMS @@ -757,6 +792,20 @@ def get_min_max (roms_data, fesom_data_lr, fesom_data_hr, x_min, x_max, y_min, y return var_min, var_max +# 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. +# Input: +# x_min, x_max, y_min, y_max = bounds on x and y (using polar coordinate +# transformation from above) for the desired region +# num_bins = number of bins to use for both x and y dimensions (try 30 or 40) +# Output: +# x_centres, y_centres = 1D arrays of length num_bins containing x and y +# coordinates of the bin centres +# roms_ubin, roms_vbin = 2D arrays of size num_bins x num_bins containing +# vector components for MetROMS at each bin +# fesom_ubin_lr, fesom_vbin_lr = same for FESOM low-res +# fesom_ubin_hr, fesom_vbin_hr = same for FESOM high-res def make_vectors (x_min, x_max, y_min, y_max, num_bins): # Set up bins (edges) @@ -847,6 +896,16 @@ def make_vectors (x_min, x_max, y_min, y_max, num_bins): return x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, fesom_ubin_hr, fesom_vbin_hr +# Plot ice shelf draft 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 ice shelf draft title, for +# use in a figure showing multiple variables def plot_draft (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 @@ -905,6 +964,24 @@ def plot_draft (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter): cbar = colorbar(img, cax=cbaxes, ticks=cbar_ticks) +# Plot ice shelf melt rate 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 +# change_points = list of size 3 containing values where the colourmap should +# transition (1) from yellow to orange, (2) from orange to red, +# (3) from red to magenta. Should not include the minimum value, +# 0, or the maximum value. This way the custom colourmap can be +# adjusted so that all melt rates are visible, particularly +# for ice shelves with strong spatial variations in melt. +# letter = 'a', 'b', 'c', etc. to add before the ice shelf melt rate title, for +# use in a figure showing multiple variables +# y0 = y-coordinate of model titles for the entire plot, assuming melt rate is +# always at the top (i.e. letter='a'). Play around between 1.15 and 1.35. def plot_melt (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, change_points, letter, y0): # Set up a grey square for FESOM to fill the background with land @@ -987,6 +1064,16 @@ def plot_melt (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, change_points text(0.5, y0, 'FESOM (high-res)', fontsize=18, horizontalalignment='center', transform=ax.transAxes) +# Plot bottom water temperature 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 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): # Set up a grey square for FESOM to fill the background with land @@ -1045,6 +1132,16 @@ def plot_bwtemp (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter): cbar = colorbar(img, cax=cbaxes, ticks=cbar_ticks) +# Plot bottom water salinity 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 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): # Set up a grey square for FESOM to fill the background with land @@ -1103,6 +1200,22 @@ def plot_bwsalt (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, letter): cbar = colorbar(img, cax=cbaxes, ticks=cbar_ticks) +# Plot vertically averaged velocity 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 +# x_centres, y_centres, roms_ubin, roms_vbin, fesom_ubin_lr, fesom_vbin_lr, +# fesom_ubin_hr, fesom_vbin_hr = output variables from the +# "make_vectors" function for the given region +# letter = 'a', 'b', 'c', etc. to add before the bottom water salinity title, +# for use in a figure showing multiple variables +# loc_string = optional string containing location title, in the case of +# velocity being zoomed into a smaller region to the other +# variables (eg "George VI Ice Shelf" for the Bellingshausen plot) 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): # Set up a grey square for FESOM to fill the background with land @@ -1169,6 +1282,24 @@ def plot_velavg (x_min, x_max, y_min, y_max, gs, cbaxes, cbar_ticks, x_centres, cbar = colorbar(img, cax=cbaxes, ticks=cbar_ticks) +# Plot zonal slices (latitude vs depth) of temperature and salinity for +# MetROMS, low-res FESOM, and high-res FESOM through the given longitude, in the +# given axes. +# Input: +# lon0 = longitude to plot (-180 to 180) +# lat_min, lat_max = bounds on latitude to plot (-90 to 90) +# depth_min, depth_max = bounds on depth to plot (negative, in metres) +# gs1, gs2 = GridSpec objects of size 1x3, in which to plot temperature and +# salinity respectively +# cbaxes1, cbaxes2 = Axes objects for locations of temperature and salinity +# colourbars +# cbar_ticks1, cbar_ticks2 = 1D arrays containing values for ticks on +# temperature and salinity colourbars +# letter1, letter2 = 'a', 'b', 'c', etc. to add before the temperature and +# salinity titles +# loc_string = optional string containing location title, in the case of +# plotting slices through a single ice shelf in a large region +# (eg "Fimbul Ice Shelf" for the Eastern Weddell plot) def plot_zonal_ts (lon0, lat_min, lat_max, depth_min, depth_max, gs1, gs2, cbaxes1, cbaxes2, cbar_ticks1, cbar_ticks2, letter1, letter2, loc_string=None): # Figure out what to write on the title about longitude diff --git a/mip_regions_1var.py b/mip_regions_1var.py index 391904d..9f2387b 100644 --- a/mip_regions_1var.py +++ b/mip_regions_1var.py @@ -17,11 +17,11 @@ # 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 -# melt rate, bottom water temperature, bottom water salinity, surface 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 surface velocity figures show absolute value -# (speed) in colours, and direction with vector arrows. +# 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. def mip_regions_1var (): # Path to ROMS grid file @@ -51,7 +51,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 = ['draft', 'melt', 'temp', 'salt', 'vel'] + var_names = ['vel'] #['draft', 'melt', 'temp', 'salt', 'vel'] # Constants sec_per_year = 365*24*3600 deg2rad = pi/180.0 @@ -73,6 +73,7 @@ def mip_regions_1var (): id = Dataset(roms_grid, 'r') roms_lon = id.variables['lon_rho'][:,:] roms_lat = id.variables['lat_rho'][:,:] + roms_h = id.variables['h'][:,:] roms_mask = id.variables['mask_rho'][:,:] roms_zice = id.variables['zice'][:,:] roms_angle = id.variables['angle'][:,:] @@ -118,37 +119,48 @@ def mip_regions_1var (): # Bottom layer roms_data = id.variables[var][0,0,:,:] elif var == 'vel': - # Read surface u and v - u_tmp = id.variables['u'][0,-1,:,:] - v_tmp = id.variables['v'][0,-1,:,:] - # Extend into land mask before interpolation to rho-grid so - # the land mask doesn't change in the final plot - num_lat_u = size(u_tmp,0) - num_lon_u = size(u_tmp,1) - for j in range(1,num_lat_u-1): - for i in range(1,num_lon_u-1): - # Check for masked points - if u_tmp[j,i] is ma.masked: - # Look at 4 neighbours - neighbours = ma.array([u_tmp[j-1,i], u_tmp[j,i-1], u_tmp[j+1,i], u_tmp[j,i+1]]) - # Find how many of them are unmasked - num_unmasked = MaskedArray.count(neighbours) - if num_unmasked > 0: - # There is at least one unmasked neighbour; - # set u_tmp to their average - u_tmp[j,i] = sum(neighbours)/num_unmasked - # Repeat for v - num_lat_v = size(v_tmp,0) - num_lon_v = size(v_tmp,1) - for j in range(1,num_lat_v-1): - for i in range(1,num_lon_v-1): - if v_tmp[j,i] is ma.masked: - neighbours = ma.array([v_tmp[j-1,i], v_tmp[j,i-1], v_tmp[j+1,i], v_tmp[j,i+1]]) - num_unmasked = MaskedArray.count(neighbours) - if num_unmasked > 0: - v_tmp[j,i] = sum(neighbours)/num_unmasked - # Interpolate to rho grid and rotate - u_rho, v_rho = rotate_vector_roms(u_tmp, v_tmp, roms_angle) + # Read full 3D u and v + u_3d_tmp = id.variables['u'][0,:,:,:] + v_3d_tmp = id.variables['v'][0,:,:,:] + # Get integrands on 3D grid; we only care about dz + dx, dy, dz, z = cartesian_grid_3d(roms_lon, roms_lat, roms_h, roms_zice, theta_s, theta_b, hc, N) + # Unrotate each vertical level + u_3d = ma.empty(shape(dz)) + v_3d = ma.empty(shape(dz)) + num_lat_u = size(u_3d_tmp,1) + num_lon_u = size(u_3d_tmp,2) + num_lat_v = size(v_3d_tmp,1) + num_lon_v = size(v_3d_tmp,2) + for k in range(N): + # Extend into land mask before interpolation to rho-grid so + # the land mask doesn't change in the final plot + for j in range(1,num_lat_u-1): + for i in range(1,num_lon_u-1): + # Check for masked points + if u_3d_tmp[k,j,i] is ma.masked: + # Look at 4 neighbours + neighbours = ma.array([u_3d_tmp[k,j-1,i], u_3d_tmp[k,j,i-1], u_3d_tmp[k,j+1,i], u_3d_tmp[k,j,i+1]]) + # Find how many of them are unmasked + num_unmasked = MaskedArray.count(neighbours) + if num_unmasked > 0: + # There is at least one unmasked neighbour; + # set u_3d_tmp to their average + u_3d_tmp[k,j,i] = sum(neighbours)/num_unmasked + # Repeat for v + for j in range(1,num_lat_v-1): + for i in range(1,num_lon_v-1): + if v_3d_tmp[k,j,i] is ma.masked: + neighbours = ma.array([v_3d_tmp[k,j-1,i], v_3d_tmp[k,j,i-1], v_3d_tmp[k,j+1,i], v_3d_tmp[k,j,i+1]]) + num_unmasked = MaskedArray.count(neighbours) + if num_unmasked > 0: + v_3d_tmp[k,j,i] = sum(neighbours)/num_unmasked + # Interpolate to rho grid and rotate + u_k, v_k = rotate_vector_roms(u_3d_tmp[k,:,:], v_3d_tmp[k,:,:], roms_angle) + u_3d[k,:,:] = u_k + v_3d[k,:,:] = v_k + # Vertically average u and v + u_rho = sum(u_3d*dz, axis=0)/sum(dz, axis=0) + v_rho = sum(v_3d*dz, axis=0)/sum(dz, axis=0) # Get speed roms_data = sqrt(u_rho**2 + v_rho**2) id.close() @@ -185,35 +197,70 @@ def mip_regions_1var (): f.close() # Save the number of 2D nodes fesom_n2d_lr = len(fesom_cavity_lr) - # Read rotated lat and lon for each node + # Read rotated lat and lon for each node, also depth f = open(fesom_mesh_path_lr + 'nod3d.out', 'r') f.readline() rlon_lr = [] rlat_lr = [] + node_depth_lr = [] for line in f: tmp = line.split() lon_tmp = float(tmp[1]) lat_tmp = float(tmp[2]) + node_depth_tmp = -1*float(tmp[3]) if lon_tmp < -180: lon_tmp += 360 elif lon_tmp > 180: lon_tmp -= 360 rlon_lr.append(lon_tmp) rlat_lr.append(lat_tmp) + node_depth_lr.append(node_depth_tmp) f.close() # For lat and lon, only care about the 2D nodes (the first # fesom_n2d indices) rlon_lr = array(rlon_lr[0:fesom_n2d_lr]) rlat_lr = array(rlat_lr[0:fesom_n2d_lr]) + node_depth_lr = array(node_depth_lr) # Unrotate longitude fesom_lon_lr, fesom_lat_lr = unrotate_grid(rlon_lr, rlat_lr) # Calculate polar coordinates of each node fesom_x_lr = -(fesom_lat_lr+90)*cos(fesom_lon_lr*deg2rad+pi/2) fesom_y_lr = (fesom_lat_lr+90)*sin(fesom_lon_lr*deg2rad+pi/2) + # Read lists of which nodes are directly below which + f = open(fesom_mesh_path_lr + 'aux3d.out', 'r') + max_num_layers_lr = int(f.readline()) + node_columns_lr = zeros([fesom_n2d_lr, max_num_layers_lr]) + for n in range(fesom_n2d_lr): + for k in range(max_num_layers_lr): + node_columns_lr[n,k] = int(f.readline()) + node_columns_lr = node_columns_lr.astype(int) + f.close() # Now we can actually read the data - # Only care about the first fesom_n2d nodes (surface) - node_ur_lr = id.variables['u'][0,0:fesom_n2d_lr] - node_vr_lr = id.variables['v'][0,0:fesom_n2d_lr] + # Read full 3D field for both u and v + node_ur_3d_lr = id.variables['u'][0,:] + node_vr_3d_lr = id.variables['v'][0,:] + # Vertically average + node_ur_lr = zeros(fesom_n2d_lr) + node_vr_lr = zeros(fesom_n2d_lr) + for n in range(fesom_n2d_lr): + # Integrate udz, vdz, and dz over this water column + udz_col = 0 + vdz_col = 0 + dz_col = 0 + for k in range(max_num_layers_lr-1): + if node_columns_lr[n,k+1] == -999: + # Reached the bottom + break + # Trapezoidal rule + top_id = node_columns_lr[n,k] + bot_id = node_columns_lr[n,k+1] + dz_tmp = node_depth_lr[bot_id-1] - node_depth_lr[top_id-1] + udz_col += 0.5*(node_ur_3d_lr[top_id-1]+node_ur_3d_lr[bot_id-1])*dz_tmp + vdz_col += 0.5*(node_vr_3d_lr[top_id-1]+node_vr_3d_lr[bot_id-1])*dz_tmp + dz_col += dz_tmp + # Convert from integrals to averages + node_ur_lr[n] = udz_col/dz_col + node_vr_lr[n] = vdz_col/dz_col # Unrotate node_u_lr, node_v_lr = unrotate_vector(rlon_lr, rlat_lr, node_ur_lr, node_vr_lr) # Calculate speed @@ -263,24 +310,53 @@ def mip_regions_1var (): f.readline() rlon_hr = [] rlat_hr = [] + node_depth_hr = [] for line in f: tmp = line.split() lon_tmp = float(tmp[1]) lat_tmp = float(tmp[2]) + node_depth_tmp = -1*float(tmp[3]) if lon_tmp < -180: lon_tmp += 360 elif lon_tmp > 180: lon_tmp -= 360 rlon_hr.append(lon_tmp) rlat_hr.append(lat_tmp) + node_depth_hr.append(node_depth_tmp) f.close() rlon_hr = array(rlon_hr[0:fesom_n2d_hr]) rlat_hr = array(rlat_hr[0:fesom_n2d_hr]) + node_depth_hr = array(node_depth_hr) fesom_lon_hr, fesom_lat_hr = unrotate_grid(rlon_hr, rlat_hr) fesom_x_hr = -(fesom_lat_hr+90)*cos(fesom_lon_hr*deg2rad+pi/2) fesom_y_hr = (fesom_lat_hr+90)*sin(fesom_lon_hr*deg2rad+pi/2) - node_ur_hr = id.variables['u'][0,0:fesom_n2d_hr] - node_vr_hr = id.variables['v'][0,0:fesom_n2d_hr] + f = open(fesom_mesh_path_hr + 'aux3d.out', 'r') + max_num_layers_hr = int(f.readline()) + node_columns_hr = zeros([fesom_n2d_hr, max_num_layers_hr]) + for n in range(fesom_n2d_hr): + for k in range(max_num_layers_hr): + node_columns_hr[n,k] = int(f.readline()) + node_columns_hr = node_columns_hr.astype(int) + f.close() + node_ur_3d_hr = id.variables['u'][0,:] + node_vr_3d_hr = id.variables['v'][0,:] + node_ur_hr = zeros(fesom_n2d_hr) + node_vr_hr = zeros(fesom_n2d_hr) + for n in range(fesom_n2d_hr): + udz_col = 0 + vdz_col = 0 + dz_col = 0 + for k in range(max_num_layers_hr-1): + if node_columns_hr[n,k+1] == -999: + break + top_id = node_columns_hr[n,k] + bot_id = node_columns_hr[n,k+1] + dz_tmp = node_depth_hr[bot_id-1] - node_depth_hr[top_id-1] + udz_col += 0.5*(node_ur_3d_hr[top_id-1]+node_ur_3d_hr[bot_id-1])*dz_tmp + vdz_col += 0.5*(node_vr_3d_hr[top_id-1]+node_vr_3d_hr[bot_id-1])*dz_tmp + dz_col += dz_tmp + node_ur_hr[n] = udz_col/dz_col + node_vr_hr[n] = vdz_col/dz_col node_u_hr, node_v_hr = unrotate_vector(rlon_hr, rlat_hr, node_ur_hr, node_vr_hr) node_data_hr = sqrt(node_u_hr**2 + node_v_hr**2) id.close() @@ -294,226 +370,226 @@ def mip_regions_1var (): 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 = ' surface 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 == '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_zonal_cavity_ts.py b/mip_zonal_cavity_ts.py index c0c1b35..dfb1a1f 100644 --- a/mip_zonal_cavity_ts.py +++ b/mip_zonal_cavity_ts.py @@ -11,15 +11,17 @@ from fesom_grid import * from fesom_sidegrid import * -# Create one 2x2 plot for each major ice shelf: zonal slices of temperature -# (top) and salinity (bottom), comparing MetROMS (left) and FESOM (right). -# Longitudes to slice through, and latitude bounds, are pre-determined. +# Create one 3x2 plot for each major ice shelf: zonal slices of temperature +# (top) and salinity (bottom), comparing MetROMS, low-res FESOM, and high-res +# FESOM. Longitudes to slice through, and latitude bounds, are pre-determined. # Input: # roms_grid = path to ROMS grid file # roms_file = path to ROMS time-averaged file containing 3D temp and salt -# fesom_mesh_path = path to FESOM mesh directory -# fesom_file = path to FESOM time-averaged file containing 3D temp and salt -def mip_zonal_cavity_ts (roms_grid, roms_file, fesom_mesh_path, fesom_file): +# fesom_mesh_path_lr, fesom_mesh_path_hr = paths to FESOM mesh directories for +# the low-res and high-res simulations +# fesom_file_lr, fesom_file_hr = paths to FESOM time-averaged files containing +# 3D temp and salt for the low-res and high-res simulations +def mip_zonal_cavity_ts (roms_grid, roms_file, fesom_mesh_path_lr, fesom_file_lr, fesom_mesh_path_hr, fesom_file_hr): # Name of each ice shelf shelf_names = ['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'] @@ -53,13 +55,20 @@ def mip_zonal_cavity_ts (roms_grid, roms_file, fesom_mesh_path, fesom_file): roms_salt_3d = id.variables['salt'][0,:,:,:] id.close() - print 'Setting up FESOM' + print 'Setting up low-res FESOM' # Build the regular FESOM grid - elm2D = fesom_grid(fesom_mesh_path) + elm2D_lr = fesom_grid(fesom_mesh_path_lr) # Read temperature and salinity at every node - id = Dataset(fesom_file, 'r') - fesom_temp_nodes = id.variables['temp'][0,:] - fesom_salt_nodes = id.variables['salt'][0,:] + id = Dataset(fesom_file_lr, 'r') + fesom_temp_nodes_lr = id.variables['temp'][0,:] + fesom_salt_nodes_lr = id.variables['salt'][0,:] + id.close() + + print 'Setting up high-res FESOM' + elm2D_hr = fesom_grid(fesom_mesh_path_hr) + id = Dataset(fesom_file_hr, 'r') + fesom_temp_nodes_hr = id.variables['temp'][0,:] + fesom_salt_nodes_hr = id.variables['salt'][0,:] id.close() # Loop over ice shelves @@ -85,78 +94,114 @@ def mip_zonal_cavity_ts (roms_grid, roms_file, fesom_mesh_path, fesom_file): # Round down to nearest 50 metres depth_min = floor(depth_min_tmp/50)*50 - # FESOM + # FESOM low-res # Build arrays of SideElements making up zonal slices - selements_temp = fesom_sidegrid(elm2D, fesom_temp_nodes, lon0[index], lat_max[index]) - selements_salt = fesom_sidegrid(elm2D, fesom_salt_nodes, lon0[index], lat_max[index]) + selements_temp_lr = fesom_sidegrid(elm2D_lr, fesom_temp_nodes_lr, lon0[index], lat_max[index]) + selements_salt_lr = fesom_sidegrid(elm2D_lr, fesom_salt_nodes_lr, lon0[index], lat_max[index]) # Build array of quadrilateral patches for the plots, and data values # corresponding to each SideElement - patches = [] - fesom_temp = [] - for selm in selements_temp: + patches_lr = [] + fesom_temp_lr = [] + for selm in selements_temp_lr: # Make patch coord = transpose(vstack((selm.y, selm.z))) - patches.append(Polygon(coord, True, linewidth=0.)) + patches_lr.append(Polygon(coord, True, linewidth=0.)) # Save data value - fesom_temp.append(selm.var) - fesom_temp = array(fesom_temp) + fesom_temp_lr.append(selm.var) + fesom_temp_lr = array(fesom_temp_lr) # Salinity has same patches but different values - fesom_salt = [] - for selm in selements_salt: - fesom_salt.append(selm.var) - fesom_salt = array(fesom_salt) + fesom_salt_lr = [] + for selm in selements_salt_lr: + fesom_salt_lr.append(selm.var) + fesom_salt_lr = array(fesom_salt_lr) + + # FESOM high-res + selements_temp_hr = fesom_sidegrid(elm2D_hr, fesom_temp_nodes_hr, lon0[index], lat_max[index]) + selements_salt_hr = fesom_sidegrid(elm2D_hr, fesom_salt_nodes_hr, lon0[index], lat_max[index]) + patches_hr = [] + fesom_temp_hr = [] + for selm in selements_temp_hr: + coord = transpose(vstack((selm.y, selm.z))) + patches_hr.append(Polygon(coord, True, linewidth=0.)) + fesom_temp_hr.append(selm.var) + fesom_temp_hr = array(fesom_temp_hr) + fesom_salt_hr = [] + for selm in selements_salt_hr: + fesom_salt_hr.append(selm.var) + fesom_salt_hr = array(fesom_salt_hr) # Find bounds on each variable - temp_min = min(amin(roms_temp[flag]), amin(fesom_temp)) - temp_max = max(amax(roms_temp[flag]), amax(fesom_temp)) - salt_min = min(amin(roms_salt[flag]), amin(fesom_salt)) - salt_max = max(amax(roms_salt[flag]), amax(fesom_salt)) + temp_min = amin(array([amin(roms_temp[flag]), amin(fesom_temp_lr), amin(fesom_temp_hr)])) + temp_max = amax(array([amax(roms_temp[flag]), amax(fesom_temp_lr), amax(fesom_temp_hr)])) + salt_min = amin(array([amin(roms_salt[flag]), amin(fesom_salt_lr), amin(fesom_salt_hr)])) + salt_max = amax(array([amax(roms_salt[flag]), amax(fesom_salt_lr), amax(fesom_salt_hr)])) # Plot - fig = figure(figsize=(18,12)) + fig = figure(figsize=(24,12)) # MetROMS temperature - ax = fig.add_subplot(2, 2, 1) + ax = fig.add_subplot(2, 3, 1) pcolor(roms_lat, roms_z, roms_temp, vmin=temp_min, vmax=temp_max, cmap='jet') title(r'MetROMS temperature ($^{\circ}$C)', fontsize=20) ylabel('Depth (m)', fontsize=16) xlim([lat_min[index], lat_max[index]]) ylim([depth_min, 0]) - # FESOM temperature - ax = fig.add_subplot(2, 2, 2) - img1 = PatchCollection(patches, cmap='jet') - img1.set_array(fesom_temp) - img1.set_edgecolor('face') - img1.set_clim(vmin=temp_min, vmax=temp_max) - ax.add_collection(img1) - title(r'FESOM temperature ($^{\circ}$C)', fontsize=20) + # FESOM low-res temperature + ax = fig.add_subplot(2, 3, 2) + img = PatchCollection(patches_lr, cmap='jet') + img.set_array(fesom_temp_lr) + img.set_edgecolor('face') + img.set_clim(vmin=temp_min, vmax=temp_max) + ax.add_collection(img) + title(r'FESOM (low-res) temperature ($^{\circ}$C)', fontsize=20) + xlim([lat_min[index], lat_max[index]]) + ylim([depth_min, 0]) + # FESOM high-res temperature + ax = fig.add_subplot(2, 3, 3) + img = PatchCollection(patches_hr, cmap='jet') + img.set_array(fesom_temp_hr) + img.set_edgecolor('face') + img.set_clim(vmin=temp_min, vmax=temp_max) + ax.add_collection(img) + title(r'FESOM (high-res) temperature ($^{\circ}$C)', fontsize=20) xlim([lat_min[index], lat_max[index]]) ylim([depth_min, 0]) # Add colorbar for temperature - cbaxes1 = fig.add_axes([0.92, 0.575, 0.01, 0.3]) - cbar1 = colorbar(img1, cax=cbaxes1) - cbar1.ax.tick_params(labelsize=16) + cbaxes = fig.add_axes([0.92, 0.575, 0.01, 0.3]) + cbar = colorbar(img, cax=cbaxes) + cbar.ax.tick_params(labelsize=16) # MetROMS salinity - ax = fig.add_subplot(2, 2, 3) + ax = fig.add_subplot(2, 3, 4) pcolor(roms_lat, roms_z, roms_salt, vmin=salt_min, vmax=salt_max, cmap='jet') title('MetROMS salinity (psu)', fontsize=20) xlabel('Latitude', fontsize=16) ylabel('Depth (m)', fontsize=16) xlim([lat_min[index], lat_max[index]]) ylim([depth_min, 0]) - # FESOM salinity - ax = fig.add_subplot(2, 2, 4) - img2 = PatchCollection(patches, cmap='jet') - img2.set_array(fesom_salt) - img2.set_edgecolor('face') - img2.set_clim(vmin=salt_min, vmax=salt_max) - ax.add_collection(img2) - title(r'FESOM salinity (psu)', fontsize=20) + # FESOM low-res salinity + ax = fig.add_subplot(2, 3, 5) + img = PatchCollection(patches_lr, cmap='jet') + img.set_array(fesom_salt_lr) + img.set_edgecolor('face') + img.set_clim(vmin=salt_min, vmax=salt_max) + ax.add_collection(img) + title(r'FESOM (low-res) salinity (psu)', fontsize=20) + xlabel('Latitude', fontsize=16) + xlim([lat_min[index], lat_max[index]]) + ylim([depth_min, 0]) + # FESOM high-res salinity + ax = fig.add_subplot(2, 3, 6) + img = PatchCollection(patches_hr, cmap='jet') + img.set_array(fesom_salt_hr) + img.set_edgecolor('face') + img.set_clim(vmin=salt_min, vmax=salt_max) + ax.add_collection(img) + title(r'FESOM (high-res) salinity (psu)', fontsize=20) xlabel('Latitude', fontsize=16) xlim([lat_min[index], lat_max[index]]) ylim([depth_min, 0]) # Add colorbar for salinity - cbaxes2 = fig.add_axes([0.92, 0.125, 0.01, 0.3]) - cbar2 = colorbar(img2, cax=cbaxes2) - cbar2.ax.tick_params(labelsize=16) + cbaxes = fig.add_axes([0.92, 0.125, 0.01, 0.3]) + cbar = colorbar(img, cax=cbaxes) + cbar.ax.tick_params(labelsize=16) # Main title suptitle(shelf_names[index] + lon_string, fontsize=28) #fig.show() @@ -168,6 +213,8 @@ def mip_zonal_cavity_ts (roms_grid, roms_file, fesom_mesh_path, fesom_file): roms_grid = raw_input("Path to ROMS grid file: ") roms_file = raw_input("Path to ROMS time-averaged file containing 3D temp and salt: ") - fesom_mesh_path = raw_input("Path to FESOM mesh directory: ") - fesom_file = raw_input("Path to FESOM time-averaged file containing 3D temp and salt: ") - mip_zonal_cavity_ts(roms_grid, roms_file, fesom_mesh_path, fesom_file) + fesom_mesh_path_lr = raw_input("Path to FESOM low-res mesh directory: ") + fesom_file_lr = raw_input("Path to FESOM low-res time-averaged file containing 3D temp and salt: ") + fesom_mesh_path_hr = raw_input("Path to FESOM high-res mesh directory: ") + fesom_file_hr = raw_input("Path to FESOM high-res time-averaged file containing 3D temp and salt: ") + mip_zonal_cavity_ts(roms_grid, roms_file, fesom_mesh_path_lr, fesom_file_lr, fesom_mesh_path_hr, fesom_file_hr)