Final scripts for first draft of thesis

bugs_acc_fig.py

@@ -1,22 +1,50 @@
 from netCDF4 import Dataset, num2date
 from numpy import *
 from matplotlib.pyplot import *
+from matplotlib.colors import ListedColormap, LinearSegmentedColormap
 from rotate_vector_roms import *
 
+# Used ocean_his_0102.nc, timestep 1
 def bugs_acc_fig (grid_path, laplacian_file, biharmonic_file, tstep):
 
     # Month names for title
     month_names = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October', 'November', 'December']
     # Degrees to radians conversion factor
     deg2rad = pi/180.0
+    # Bounds on plot (in polar coordinate transformation)
+    x_min = -40
+    x_max = 0
+    y_min = 0
+    y_max = 40
+    # Minimum speed to plot (mask out values below)
+    threshold = 0.12
+    # Maximum speed to plot
+    bound = 1
 
     # Read angle and grid
     id = Dataset(grid_path, 'r')
     angle = id.variables['angle'][:-15,:]
     lon = id.variables['lon_rho'][:-15,:-1]
     lat = id.variables['lat_rho'][:-15,:-1]
+    mask = id.variables['mask_rho'][:-15,:-1]
     id.close()
 
+    # Set up grey shading of land
+    mask = ma.masked_where(mask==1, mask)
+    grey_cmap = ListedColormap([(0.6, 0.6, 0.6)])
+    x_reg, y_reg = meshgrid(linspace(x_min, x_max, num=1000), linspace(y_min, y_max, num=1000))
+    land_circle = zeros(shape(x_reg))
+    lon_c = 50
+    lat_c = -83
+    radius = 10.1
+    x_c = -(lat_c+90)*cos(lon_c*deg2rad+pi/2)
+    y_c = (lat_c+90)*sin(lon_c*deg2rad+pi/2)
+    land_circle = ma.masked_where(sqrt((x_reg-x_c)**2 + (y_reg-y_c)**2) > radius, land_circle)
+    # Truncate colourmap
+    min_colour = threshold/bound
+    max_colour = 1
+    trunc_cmap = truncate_colormap(get_cmap('jet'), min_colour, max_colour)
+
     # Read surface velocity
     id = Dataset(laplacian_file, 'r')
     ur_lap = id.variables['u'][tstep-1,-1,:-15,:]
@@ -42,35 +70,50 @@ def bugs_acc_fig (grid_path, laplacian_file, biharmonic_file, tstep):
     # Calculate x and y coordinates for plotting circumpolar projection
     x = -(lat+90)*cos(lon*deg2rad+pi/2)
     y  = (lat+90)*sin(lon*deg2rad+pi/2)
+    speed_laplacian = ma.masked_where(speed_laplacian < threshold, speed_laplacian)
+    speed_biharmonic = ma.masked_where(speed_biharmonic < threshold, speed_biharmonic)
 
     # Make the plot
     fig = figure(figsize=(20,10))
     # Laplacian
     ax = fig.add_subplot(1,2,1, aspect='equal')
-    pcolor(x, y, speed_laplacian, vmin=0, vmax=1.25, cmap='cool')
+    # First shade land in grey
+    pcolor(x, y, mask, cmap=grey_cmap)
+    pcolor(x_reg, y_reg, land_circle, cmap=grey_cmap)
+    pcolor(x, y, speed_laplacian, vmin=threshold, vmax=bound, cmap=trunc_cmap)
     title('Laplacian viscosity', fontsize=24)
-    xlim([-40, 0])
-    ylim([0, 40])
+    xlim([x_min, x_max])
+    ylim([y_min, y_max])
     ax.set_xticks([])
     ax.set_yticks([])
     # Biharmonic
     ax = fig.add_subplot(1,2,2, aspect='equal')
-    img = pcolor(x, y, speed_biharmonic, vmin=0, vmax=1.25, cmap='cool')
+    pcolor(x, y, mask, cmap=grey_cmap)
+    pcolor(x_reg, y_reg, land_circle, cmap=grey_cmap)
+    img = pcolor(x, y, speed_biharmonic, vmin=threshold, vmax=bound, cmap=trunc_cmap)
     title('Biharmonic viscosity', fontsize=24)
     xlim([-40, 0])
     ylim([0, 40])
     ax.set_xticks([])
     ax.set_yticks([])
     # Colourbar on the right
     cbaxes = fig.add_axes([0.94, 0.3, 0.02, 0.4])
-    cbar = colorbar(img, cax=cbaxes, extend='max', ticks=arange(0,1.25+0.25,0.25))
+    cbar = colorbar(img, cax=cbaxes, extend='max', ticks=arange(0.2,1+0.2,0.2))
     cbar.ax.tick_params(labelsize=16)
     suptitle('Surface speed (m/s): snapshot on ' + date_string, fontsize=30)
     subplots_adjust(wspace=0.05)
     fig.show()
     fig.savefig('bugs_acc.png')
 
 
+# Truncate colourmap function from https://stackoverflow.com/questions/40929467/how-to-use-and-plot-only-a-part-of-a-colorbar-in-matplotlib
+def truncate_colormap(cmap, minval=0.0, maxval=1.0, n=-1):
+    if n== -1:
+        n = cmap.N
+    new_cmap = LinearSegmentedColormap.from_list('trunc({name},{a:.2f},{b:.2f})'.format(name=cmap.name, a=minval, b=maxval), cmap(linspace(minval, maxval, n)))
+    return new_cmap
+
+
 # Command-line interface
 if __name__ == "__main__":
 

bugs_convection_slices.py

@@ -0,0 +1,160 @@
+from netCDF4 import Dataset, num2date
+from numpy import *
+from matplotlib.pyplot import *
+from calc_z import *
+from interp_lon_roms import *
+
+def bugs_convection_slices ():
+
+    # Files and timesteps to plot
+    file_beg = '/short/m68/kaa561/metroms_iceshelf/tmproms/run/bug_chapter/no_restoring/ocean_avg_0003.nc'
+    tstep_beg = 18
+    file_end = '/short/m68/kaa561/metroms_iceshelf/tmproms/run/bug_chapter/no_restoring/ocean_avg_0012.nc'
+    tstep_end = 2
+    # Longitude to interpolate to
+    lon0 = -13
+    # Deepest depth to plot
+    depth_min = -250
+    depth_ticks = arange(depth_min, 0+50, 50)
+    depth_labels = []
+    for depth in depth_ticks:
+        depth_labels.append(str(int(-depth)))
+    # Latitudes to plot
+    lat_min = -74
+    lat_max = -55
+    lat_ticks = arange(lat_min+4, lat_max+5, 5)
+    lat_labels = []
+    for lat in lat_ticks:
+        lat_labels.append(str(int(-lat)) + r'$^{\circ}$S')
+    # Bounds on colour scales
+    var_min = [-2, 33.9]
+    var_max = [2, 34.7]
+    var_tick = [1, 0.2]
+    lev1 = linspace(var_min[0], var_max[0], num=50)
+    lev2 = linspace(var_min[1], var_max[1], num=50)
+    # Grid parameters
+    theta_s = 7.0
+    theta_b = 2.0
+    hc = 250
+    N = 31
+    # Month names for titles
+    month_names = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October', 'November', 'December']
+
+    # Get longitude for the title
+    if lon0 < 0:
+        lon_string = str(int(round(-lon0))) + r'$^{\circ}$W'
+    else:
+        lon_string = str(int(round(lon0))) + r'$^{\circ}$E'
+    # Make sure we are in the range 0-360
+    if lon0 < 0:
+        lon0 += 360
+
+    # Set up figure
+    fig = figure(figsize=(18,12))
+    gs = GridSpec(2,2)
+    gs.update(left=0.13, right=0.9, bottom=0.05, top=0.9, wspace=0.05, hspace=0.28)
+
+    # Read 3D temperature, salinity, and grid variables at the beginning
+    id = Dataset(file_beg, 'r')
+    temp_3d_beg = id.variables['temp'][tstep_beg-1,:,:,:]
+    salt_3d_beg = id.variables['salt'][tstep_beg-1,:,:,:]
+    zeta_beg = id.variables['zeta'][tstep_beg-1,:,:]
+    h = id.variables['h'][:,:]
+    zice = id.variables['zice'][:,:]
+    lon_2d = id.variables['lon_rho'][:,:]
+    lat_2d = id.variables['lat_rho'][:,:]
+    # Read time axis and convert to Date objects
+    time_id = id.variables['ocean_time']
+    time_beg = num2date(time_id[tstep_beg-1], units=time_id.units, calendar=time_id.calendar.lower())
+    id.close()
+    # Get a 3D array of z-coordinates
+    z_3d, sc_r, Cs_r = calc_z(h, zice, theta_s, theta_b, hc, N, zeta_beg)
+    # Get the date for the title
+    date_string_beg = str(time_beg.day) + ' ' + month_names[time_beg.month-1] + ' ' + str(time_beg.year)
+    # Interpolate to lon0
+    temp_beg, z, lat = interp_lon_roms(temp_3d_beg, z_3d, lat_2d, lon_2d, lon0)
+    salt_beg, z, lat = interp_lon_roms(salt_3d_beg, z_3d, lat_2d, lon_2d, lon0)
+    # Plot temperature
+    ax = subplot(gs[0,0])
+    img = contourf(lat, z, temp_beg, lev1, extend='both', cmap='jet')
+    title(r'Temperature ($^{\circ}$C)', fontsize=24)
+    ax.set_xlim([lat_min, lat_max])
+    ax.set_xticks(lat_ticks)
+    ax.set_xticklabels([])
+    ax.set_ylim([depth_min, 0])
+    ax.set_yticks(depth_ticks)
+    ax.set_yticklabels(depth_labels, fontsize=16)
+    ylabel('Depth (m)', fontsize=18)
+    # Plot salinity
+    ax = subplot(gs[0,1])
+    img = contourf(lat, z, salt_beg, lev2, extend='both', cmap='jet')
+    title('Salinity (psu)', fontsize=24)
+    ax.set_xlim([lat_min, lat_max])
+    ax.set_xticks(lat_ticks)
+    ax.set_xticklabels([])
+    ax.set_ylim([depth_min, 0])
+    ax.set_yticks(depth_ticks)
+    ax.set_yticklabels([])
+    # Label the timestep
+    text(0.5, 0.95, date_string_beg + ', ' + lon_string, ha='center', transform=fig.transFigure, fontsize=28)
+
+    # Read 3D temperature, salinity, and sea surface height at the end
+    id = Dataset(file_end, 'r')
+    temp_3d_end = id.variables['temp'][tstep_end-1,:,:,:]
+    salt_3d_end = id.variables['salt'][tstep_end-1,:,:,:]
+    zeta_end = id.variables['zeta'][tstep_end-1,:,:]
+    # Read time axis and convert to Date objects
+    time_id = id.variables['ocean_time']
+    time_end = num2date(time_id[tstep_end-1], units=time_id.units, calendar=time_id.calendar.lower())
+    id.close()
+    # Get a 3D array of z-coordinates
+    z_3d, sc_r, Cs_r = calc_z(h, zice, theta_s, theta_b, hc, N, zeta_end)
+    # Get the date for the title
+    date_string_end = str(time_end.day) + ' ' + month_names[time_end.month-1] + ' ' + str(time_end.year)
+    # Interpolate to lon0
+    temp_end, z, lat = interp_lon_roms(temp_3d_end, z_3d, lat_2d, lon_2d, lon0)
+    salt_end, z, lat = interp_lon_roms(salt_3d_end, z_3d, lat_2d, lon_2d, lon0)
+    # Plot temperature
+    ax = subplot(gs[1,0])
+    img = contourf(lat, z, temp_end, lev1, extend='both', cmap='jet')
+    title(r'Temperature $^{\circ}$C)', fontsize=24)
+    ax.set_xlim([lat_min, lat_max])
+    ax.set_xticks(lat_ticks)
+    ax.set_xticklabels(lat_labels, fontsize=16)
+    xlabel('Latitude', fontsize=18)
+    ax.set_ylim([depth_min, 0])
+    ax.set_yticks(depth_ticks)
+    ax.set_yticklabels(depth_labels, fontsize=16)
+    ylabel('Depth (m)', fontsize=18)
+    # Colourbar
+    cbaxes = fig.add_axes([0.03, 0.3, 0.02, 0.4])
+    cbar = colorbar(img, cax=cbaxes, ticks=arange(var_min[0], var_max[0]+var_tick[0], var_tick[0]))
+    cbar.ax.tick_params(labelsize=16)
+    # Plot salinity
+    ax = subplot(gs[1,1])
+    img = contourf(lat, z, salt_end, lev2, extend='both', cmap='jet')
+    title('Salinity (psu)', fontsize=24)
+    ax.set_xlim([lat_min, lat_max])
+    ax.set_xticks(lat_ticks)
+    ax.set_xticklabels(lat_labels, fontsize=16)
+    xlabel('Latitude', fontsize=18)
+    ax.set_ylim([depth_min, 0])
+    ax.set_yticks(depth_ticks)
+    ax.set_yticklabels([])
+    # Colourbar
+    cbaxes = fig.add_axes([0.93, 0.3, 0.02, 0.4])
+    cbar = colorbar(img, cax=cbaxes, ticks=arange(var_min[1], var_max[1]+var_tick[1], var_tick[1]))
+    cbar.ax.tick_params(labelsize=16)
+    # Label the timestep
+    text(0.5, 0.47, date_string_end + ', ' + lon_string, ha='center', transform=fig.transFigure, fontsize=28)
+    fig.show()
+    fig.savefig('bugs_convection_slices.png')
+
+
+# Command-line interface
+if __name__ == "__main__":
+
+    bugs_convection_slices()
+
+
+

bugs_dpt_fig.py

@@ -47,19 +47,16 @@ def bugs_dpt_fig (dpt_log_laplacian, dpt_log_biharmonic):
     # Make new time array
     time_annual = arange(len(dpt_laplacian_avg)) + year_start
     # Plot
-    fig, ax = subplots(figsize=(10,6))
+    fig, ax = subplots(figsize=(8,6))
     ax.plot(time_annual, dpt_laplacian_avg, label='Laplacian', color='green', linewidth=2)
     ax.plot(time_annual, dpt_biharmonic_avg, label='Biharmonic', color='blue', linewidth=2)
     title('Drake Passage Transport (annually averaged)', fontsize=18)
     xlabel('Year', fontsize=14)
     ylabel('Sv', fontsize=14)
     xlim([time_annual[0], time_annual[-1]])
     grid(True)
-    # Move plot over to make room for legend
-    box = ax.get_position()
-    ax.set_position([box.x0, box.y0, box.width*0.8, box.height])
     # Make legend
-    ax.legend(loc='center left', bbox_to_anchor=(1,0.5))
+    ax.legend(loc='center right')
     fig.show()
     fig.savefig('bugs_dpt.png')
 

bugs_mld_polynya.py

@@ -0,0 +1,99 @@
+from netCDF4 import Dataset, num2date
+from numpy import *
+from matplotlib.pyplot import *
+
+def bugs_mld_polynya ():
+
+    # Files and timesteps to read
+    file_ocn = '/short/m68/kaa561/metroms_iceshelf/tmproms/run/bug_chapter/no_restoring/ocean_avg_0012.nc'
+    tstep_ocn = 2
+    file_ice = '/short/m68/kaa561/metroms_iceshelf/tmproms/run/bug_chapter/no_restoring/cice/rundir/history/iceh.1994-09-26.nc'
+    tstep_ice = 1
+    # Degrees to radians conversion factor
+    deg2rad = pi/180
+    # Month names for titles
+    month_names = ['January', 'February', 'March', 'April', 'May', 'June', 'July', 'August', 'September', 'October', 'November', 'December']
+    # Maximum latitude to plot
+    nbdry = -52+90    
+
+    # Set up figure
+    fig = figure(figsize=(16,8))
+    gs = GridSpec(1,2)
+    gs.update(left=0.1, right=0.9, bottom=0.05, top=0.9, wspace=0.05)
+
+    # Read mixed layer depth (actually surface boundary layer depth)
+    # and ROMS grid
+    id = Dataset(file_ocn, 'r')
+    lon = id.variables['lon_rho'][:,:-1]
+    lat = id.variables['lat_rho'][:,:-1]
+    zice = id.variables['zice'][:,:-1]
+    hsbl = id.variables['Hsbl'][tstep_ocn-1,:,:-1]
+    # Also read time axis and convert to Date objects
+    time_id = id.variables['ocean_time']
+    time = num2date(time_id[tstep_ocn-1], units=time_id.units, calendar=time_id.calendar.lower())
+    id.close()
+    # Mask out the ice shelves and change the sign
+    mld = ma.masked_where(zice!=0, -hsbl)
+    # Polar coordinate transformation
+    x = -(lat+90)*cos(lon*deg2rad+pi/2)
+    y = (lat+90)*sin(lon*deg2rad+pi/2)
+    # Get the date for the title
+    date_string = str(time.day) + ' ' + month_names[time.month-1] + ' ' + str(time.year)
+    # Plot
+    ax = subplot(gs[0,0], aspect='equal')
+    lev = linspace(0, amax(mld), num=50)
+    contourf(x, y, mld, lev, cmap='jet')
+    xlim([-nbdry, nbdry])
+    ylim([-nbdry, nbdry])
+    ax.set_xticks([])
+    ax.set_yticks([])
+    title('a) Surface boundary layer depth (m)', fontsize=24)
+    # Colourbar
+    cbaxes = fig.add_axes([0.03, 0.3, 0.02, 0.4])
+    cbar = colorbar(cax=cbaxes, ticks=arange(0,2000+500,500))
+    cbar.ax.tick_params(labelsize=16)
+
+    # Read sea ice concentration and CICE grid
+    id = Dataset(file_ice, 'r')
+    lon_tmp = id.variables['TLON'][:,:]
+    lat_tmp = id.variables['TLAT'][:,:]
+    aice_tmp = id.variables['aice'][tstep_ice-1,:,:]
+    id.close()
+    # Wrap the periodic boundray by 1 cell
+    lon = ma.empty([size(lon_tmp,0), size(lon_tmp,1)+1])
+    lat = ma.empty([size(lat_tmp,0), size(lat_tmp,1)+1])
+    aice = ma.empty([size(aice_tmp,0), size(aice_tmp,1)+1])
+    lon[:,:-1] = lon_tmp
+    lon[:,-1] = lon_tmp[:,0]
+    lat[:,:-1] = lat_tmp
+    lat[:,-1] = lat_tmp[:,0]
+    aice[:,:-1] = aice_tmp
+    aice[:,-1] = aice_tmp[:,0]
+    # Convert to spherical coordinates
+    x = -(lat+90)*cos(lon*deg2rad+pi/2)
+    y = (lat+90)*sin(lon*deg2rad+pi/2)
+    # Plot
+    ax = subplot(gs[0,1], aspect='equal')
+    lev = linspace(0, 1, num=50)
+    contourf(x, y, aice, lev, cmap='jet')
+    xlim([-nbdry, nbdry])
+    ylim([-nbdry, nbdry])
+    ax.set_xticks([])
+    ax.set_yticks([])
+    title('b) Sea ice concentration', fontsize=24)
+    # Colourbar
+    cbaxes = fig.add_axes([0.92, 0.3, 0.02, 0.4])
+    cbar = colorbar(cax=cbaxes, ticks=arange(0,1+0.25,0.25))
+    cbar.ax.tick_params(labelsize=16)
+    # Label the timestep
+    text(0.5, 0.94, date_string, ha='center', transform=fig.transFigure, fontsize=28)
+    fig.show()
+    fig.savefig('bugs_mld_polynya.png')
+
+
+# Command-line interface
+if __name__ == "__main__":
+
+    bugs_mld_polynya()
+
+