diff --git a/circumpolar_plot.py b/circumpolar_plot.py index 883560c..438e668 100644 --- a/circumpolar_plot.py +++ b/circumpolar_plot.py @@ -89,7 +89,9 @@ def circumpolar_plot (file_path, var_name, tstep, depth_key, depth, depth_bounds u_data = id.variables[var_name.replace('v','u')][tstep-1,k,:-15,:] u_data_lonlat, v_data_lonlat = rotate_vector_roms(u_data, v_data, angle) data_full[k,:,:] = v_data_lonlat - + + id.close() + id = Dataset(grid_path, 'r') # Read grid variables h = id.variables['h'][:-15,:] zice = id.variables['zice'][:-15,:] @@ -359,7 +361,7 @@ def average_btw_depths (data_3d, z_3d, dz_3d, z_bounds): # Will need the grid file to get the angle grid_path = raw_input("Path to ROMS grid file: ") else: - grid_path = None + grid_path = raw_input("Path to ROMS grid file: ") #grid_path = None # Get index of time axis in ROMS history/averages file tstep = int(raw_input("Timestep number (starting at 1): ")) diff --git a/file_guide.txt b/file_guide.txt index ef9d923..a27a74f 100644 --- a/file_guide.txt +++ b/file_guide.txt @@ -1165,15 +1165,16 @@ mip_timeseries.py: Plot MetROMS and FESOM timeseries together, for Drake plotted). The script will create a bunch of png files. -mip_massloss_error_map.py: Make a 2x1 circumpolar Antarctic map of percentage +mip_massloss_error_map.py: Make a 3x1 circumpolar Antarctic map of percentage error in mass loss for each ice shelf, outside the - bounds given by Rignot et al. (2013), in MetROMS - (left) and FESOM (right), for the 2003-2008 average. + bounds given by Rignot et al. (2013), in MetROMS, + FESOM low-res, and FESOM high-res, for the 2003-2008 + average. Also print the values to the screen. To run: First make sure you have run timeseries_massloss.py for ROMS through to at least the end of 2008, and the equivalent timeseries_massloss.py in the "fesomtools" - repository for the FESOM simulation. Then + repository for each FESOM simulation. Then open python or ipython and type "run mip_massloss_error_map.py". The script will prompt you for the path to the ROMS @@ -1184,23 +1185,22 @@ mip_massloss_error_map.py: Make a 2x1 circumpolar Antarctic map of percentage mip_massloss_difference_map.py: Make a circumpolar Antarctic figure showing the percentage change in mass loss for each ice - shelf in the FESOM simulation with respect to - the MetROMS simulation, for the 2003-2008 - average. + shelf in the high-resolution FESOM simulation + with respect to the low-resolution FESOM + simulation, discarding the first 10 years (i.e. + 2002-2016 average). Also print the values to + the screen. To run: First make sure you have run - timeseries_massloss.py for ROMS through - to at least the end of 2008, and the - equivalent timeseries_massloss.py in - the "fesomtools" repository for the - FESOM simulation. Then open python or - ipython and type + timeseries_massloss.py in the + "fesomtools" repository for each FESOM + simulation. Then open python or ipython + and type "run mip_massloss_error_map.py". The script will prompt you for the path - to the ROMS grid file, the ROMS and - FESOM massloss logfiles, and whether - you want to save the figure (and if so, - what filename) or display it on the - screen. + to the ROMS grid file, the FESOM + massloss logfiles, and whether you want + to save the figure (and if so, what + filename) or display it on the screen. mip_cavity_fields.py: For each major ice shelf, make a 2x1 plot of the given field for MetROMS (left) and FESOM (right), zoomed into @@ -1347,23 +1347,34 @@ mip_ts_distribution.py: Make a 2x1 plot of T/S distributions south of 65S, mesh directory, and the FESOM time-averaged file. It will output a png. -mip_drift_slices.py: Make a 4x2 plot showing zonal slices through 0E of - temperature (left) and salinity (right), in MetROMS (top) - and FESOM (bottom). The top 2x2 plots show the average - over the last year of simulation, while the bottom 2x2 - plots show anomalies with respect to the first year. +mip_drift_slices.py: Make a 3x2 plot of temperature (left) and salinity (right) + through 0E. The top row is the initial conditions from + ECCO2. The middle and bottom rows are the end of the + simulation (last 5 day average) from MetROMS and FESOM + respectively. 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. Next, make - time-averaged files of temperature and salinity - for each model over the first and last year (just - use ncra). Then open python or ipython and type - "run mip_drift_slices.py". The script will prompt - you for the paths to the ROMS grid file, the ROMS - time-averaged files, the FESOM mesh directory, and - the FESOM time-averaged files. It will output a - png. + repository on your system. Also make sure the + paths to the original ECCO2 files (ecco_temp_file, + ecco_salt_file) are correct. Then open python or + ipython and type "run mip_drift_slices.py". The + script will prompt you for the paths to the ROMS + grid file, the last ocean_avg file output by ROMS, + the FESOM mesh directory, and the last oce.mean + file output by FESOM. It will output a png. + +mip_calc_massloss.py: Calculate the average mass loss for each ice shelf over + 2002-2016 for MetROMS, FESOM low-res, and FESOM high-res. + Print to the screen along with Rignot's estimates for + 2003-2008. + To run: First make sure you have run + timeseries_massloss.py for the MetROMS simulation, + as well as the equivalent timeseries_massloss.py + in the fesomtools repository for both FESOM + simulations. Then open python or ipython and type + "run mip_calc_massloss.py". The script will + prompt you for the paths to the three logfiles. ***UTILITY FUNCTIONS*** diff --git a/mip_calc_massloss.py b/mip_calc_massloss.py new file mode 100644 index 0000000..f757fe3 --- /dev/null +++ b/mip_calc_massloss.py @@ -0,0 +1,144 @@ +from numpy import * + +# Calculate the average mass loss for each ice shelf over 2002-2016 for MetROMS, +# FESOM low-res, and FESOM high-res. Print to the screen along with Rignot's +# estimates for 2003-2008. +# Input: +# roms_logfile = path to ROMS logfile from timeseries_massloss.py +# fesom_logfile_lr, fesom_logfile_hr = paths to FESOM logfiles from +# timeseries_massloss.py in the fesomtools repository, for +# low-res and high-res respectively +def mip_calc_massloss (roms_logfile, fesom_logfile_lr, fesom_logfile_hr): + + # Year simulations start + year_start = 1992 + # Years to avearge over + calc_start = 2002 + calc_end = 2016 + # Number of output steps per year in FESOM + peryear = 365/5 + # Name of each ice shelf + names = ['Total Mass Loss', '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'] + # Observed mass loss (Rignot 2013) and uncertainty for each ice shelf, in Gt/y + obs_massloss = [1325, 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 = [235, 14, 67, 40, 45, 19, 8, 7, 4, 14, 2, 3, 3, 8, 15, 10, 23, 3, 18, 2, 2, 14, 16, 34] + num_shelves = len(obs_massloss) + + # Read ROMS logfile + roms_time = [] + f = open(roms_logfile, 'r') + # Skip the first line (header for time array) + f.readline() + for line in f: + try: + roms_time.append(float(line)) + except(ValueError): + # Reached the header for the next variable + break + # Set up array for mass loss values at each ice shelf + roms_massloss_ts = empty([num_shelves, len(roms_time)]) + index = 0 + # Loop over ice shelves + while index < num_shelves: + t = 0 + for line in f: + try: + roms_massloss_ts[index, t] = float(line) + t += 1 + except(ValueError): + # Reached the header for the next ice shelf + break + index +=1 + f.close() + # Add start year to ROMS time array + roms_time = array(roms_time) + year_start + # Average between given years + t_start = nonzero(roms_time >= calc_start)[0][0] + if calc_end == 2016: + t_end = size(roms_time) + else: + t_end = nonzero(roms_time >= calc_end+1)[0][0] + roms_massloss = mean(roms_massloss_ts[:,t_start:t_end], axis=1) + + # Read FESOM timeseries + # Low-res + tmp = [] + f = open(fesom_logfile_lr, 'r') + # Skip the first line (header) + f.readline() + # Read total mass loss + num_time = 0 + for line in f: + try: + tmp.append(float(line)) + num_time += 1 + except(ValueError): + # Reached the header for the next variable + break + # Set up array for mass loss values at each ice shelf + fesom_massloss_ts_lr = empty([num_shelves, num_time]) + # Save the total values in the first index + fesom_massloss_ts_lr[0,:] = array(tmp) + # Loop over ice shelves + index = 1 + while index < num_shelves: + t = 0 + for line in f: + try: + fesom_massloss_ts_lr[index,t] = float(line) + t += 1 + except(ValueError): + # Reached the header for the next ice shelf + break + index += 1 + f.close() + # Average between given years + fesom_massloss_lr = mean(fesom_massloss_ts_lr[:,peryear*(calc_start-year_start):peryear*(calc_end+1-year_start)], axis=1) + # Repeat for high-res + tmp = [] + f = open(fesom_logfile_hr, 'r') + f.readline() + num_time = 0 + for line in f: + try: + tmp.append(float(line)) + num_time += 1 + except(ValueError): + break + fesom_massloss_ts_hr = empty([num_shelves, num_time]) + fesom_massloss_ts_hr[0,:] = array(tmp) + index = 1 + while index < num_shelves: + t = 0 + for line in f: + try: + fesom_massloss_ts_hr[index,t] = float(line) + t += 1 + except(ValueError): + break + index += 1 + f.close() + fesom_massloss_hr = mean(fesom_massloss_ts_hr[:,peryear*(calc_start-year_start):peryear*(calc_end+1-year_start)], axis=1) + + # Loop over ice shelves + for index in range(num_shelves): + print names[index] + print 'MetROMS: ' + str(roms_massloss[index]) + print 'FESOM low-res: ' + str(fesom_massloss_lr[index]) + print 'FESOM high-res: ' + str(fesom_massloss_hr[index]) + # Find the range of observations + massloss_low = obs_massloss[index] - obs_massloss_error[index] + massloss_high = obs_massloss[index] + obs_massloss_error[index] + print 'Rignot: ' + str(massloss_low) + '-' + str(massloss_high) + + +# Command-line interface +if __name__ == "__main__": + + roms_logfile = raw_input("Path to ROMS logfile from timeseries_massloss.py: ") + fesom_logfile_lr = raw_input("Path to FESOM low-res logfile from timeseries_massloss.py: ") + fesom_logfile_hr = raw_input("Path to FESOM high-res logfile from timeseries_massloss.py: ") + mip_calc_massloss(roms_logfile, fesom_logfile_lr, fesom_logfile_hr) + + + diff --git a/mip_drift_slices.py b/mip_drift_slices.py index 67ba713..99584c0 100644 --- a/mip_drift_slices.py +++ b/mip_drift_slices.py @@ -11,28 +11,26 @@ from fesom_grid import * from fesom_sidegrid import * -# Make a 4x2 plot showing zonal slices through 0E of temperature (left) and -# salinity (right), in MetROMS (top) and FESOM (bottom). The top 2x2 plots show -# the average over the last year of simulation, while the bottom 2x2 plots show -# anomalies with respect to the first year. +# Make a 3x2 plot of temperature (left) and salinity (right) through 0E. +# The top row is the initial conditions from ECCO2. The middle and bottom rows +# are the end of the simulation (last 5 day average) from MetROMS and FESOM +# respectively. # Input: # roms_grid = path to ROMS grid file -# roms_file_first, roms_file_last = paths to files containing ROMS temperature -# and salinity averaged over the first year -# and the last year respectively -# fesom_mesh_path = path to FESOM mesh directory -# fesom_file_first, fesom_file_last = paths to files containing FESOM -# temperature and salinity averaged over -# the first year and the last year -# respectively -def mip_drift_slices (roms_grid, roms_file_first, roms_file_last, fesom_mesh_path, fesom_file_first, fesom_file_last): +# roms_file = path to last ocean_avg file output by ROMS (assume 5-day averages) +# fesom_mesh_path = path to last oce.mean file output by FESOM (assume 5-day +# averages) +def mip_drift_slices (roms_grid, roms_file, fesom_mesh_path, fesom_file): - # Longitude to plot (0E) + # Paths to ECCO2 files with initial conditions for temp and salt + ecco_temp_file = '/short/m68/kaa561/metroms_iceshelf/data/originals/ECCO2/THETA.1440x720x50.199201.nc' + ecco_salt_file = '/short/m68/kaa561/metroms_iceshelf/data/originals/ECCO2/SALT.1440x720x50.199201.nc' + # Longitude to interpolate to (OE) lon0 = 0 # Bounds on plot lat_min = -73 lat_max = -30 - depth_min = -5300 + depth_min = -6000 depth_max = 0 # ROMS grid parameters theta_s = 7.0 @@ -42,17 +40,32 @@ def mip_drift_slices (roms_grid, roms_file_first, roms_file_last, fesom_mesh_pat # Bounds on colour scales for temperature and salinity temp_min = -2 temp_max = 6 - temp_anom = 2 salt_min = 33.9 salt_max = 34.9 - salt_anom = 0.5 - # Get longitude for the title + # 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' + print 'Processing ECCO2' + id = Dataset(ecco_temp_file, 'r') + # Read grid variables + ecco_lat = id.variables['LATITUDE_T'][:] + ecco_depth = -1*id.variables['DEPTH_T'][:] + if lon0 == 0: + # Hard-coded lon0 = 0E: average between the first (0.125 E) and last + # (359.875 E = -0.125 W) indices in the regular ECCO2 grid + ecco_temp = 0.5*(id.variables['THETA'][0,:,:,0] + id.variables['THETA'][0,:,:,-1]) + id.close() + id = Dataset(ecco_salt_file, 'r') + ecco_salt = 0.5*(id.variables['SALT'][0,:,:,0] + id.variables['SALT'][0,:,:,-1]) + id.close() + else: + print 'lon0 is only coded for 0E at this time' + return + print 'Processing ROMS' # Read grid variables we need id = Dataset(roms_grid, 'r') @@ -62,27 +75,18 @@ def mip_drift_slices (roms_grid, roms_file_first, roms_file_last, fesom_mesh_pat roms_zice = id.variables['zice'][:,:] id.close() # Read temperature and salinity - id = Dataset(roms_file_first, 'r') - roms_temp_first_3d = id.variables['temp'][0,:,:,:] - roms_salt_first_3d = id.variables['salt'][0,:,:,:] + id = Dataset(roms_file, 'r') + roms_temp_3d = id.variables['temp'][-1,:,:,:] + roms_salt_3d = id.variables['salt'][-1,:,:,:] id.close() - id = Dataset(roms_file_last, 'r') - roms_temp_last_3d = id.variables['temp'][0,:,:,:] - roms_salt_last_3d = id.variables['salt'][0,:,:,:] - id.close() - # Get anomalies for last year - roms_temp_anom_3d = roms_temp_last_3d - roms_temp_first_3d - roms_salt_anom_3d = roms_salt_last_3d - roms_salt_first_3d # Get a 3D array of z-coordinates; sc_r and Cs_r are unused in this script roms_z_3d, sc_r, Cs_r = calc_z(roms_h, roms_zice, theta_s, theta_b, hc, N) # Make sure we are in the range 0-360 if lon0 < 0: lon0 += 360 # Interpolate to lon0 - roms_temp, roms_z, roms_lat = interp_lon_roms(roms_temp_last_3d, roms_z_3d, roms_lat_2d, roms_lon_2d, lon0) - roms_salt, roms_z, roms_lat = interp_lon_roms(roms_salt_last_3d, roms_z_3d, roms_lat_2d, roms_lon_2d, lon0) - roms_temp_anom, roms_z, roms_lat = interp_lon_roms(roms_temp_anom_3d, roms_z_3d, roms_lat_2d, roms_lon_2d, lon0) - roms_salt_anom, roms_z, roms_lat = interp_lon_roms(roms_salt_anom_3d, roms_z_3d, roms_lat_2d, roms_lon_2d, lon0) + roms_temp, roms_z, roms_lat = interp_lon_roms(roms_temp_3d, roms_z_3d, roms_lat_2d, roms_lon_2d, lon0) + roms_salt, roms_z, roms_lat = interp_lon_roms(roms_salt_3d, roms_z_3d, roms_lat_2d, roms_lon_2d, lon0) # Switch back to range -180-180 if lon0 > 180: lon0 -= 360 @@ -91,22 +95,13 @@ def mip_drift_slices (roms_grid, roms_file_first, roms_file_last, fesom_mesh_pat # Build regular elements elements = fesom_grid(fesom_mesh_path) # Read temperature and salinity - id = Dataset(fesom_file_first, 'r') - fesom_temp_first_nodes = id.variables['temp'][0,:] - fesom_salt_first_nodes = id.variables['salt'][0,:] - id.close() - id = Dataset(fesom_file_last, 'r') - fesom_temp_last_nodes = id.variables['temp'][0,:] - fesom_salt_last_nodes = id.variables['salt'][0,:] + id = Dataset(fesom_file, 'r') + fesom_temp_nodes = id.variables['temp'][-1,:] + fesom_salt_nodes = id.variables['salt'][-1,:] id.close() - # Get anomalies for last year - fesom_temp_anom_nodes = fesom_temp_last_nodes - fesom_temp_first_nodes - fesom_salt_anom_nodes = fesom_salt_last_nodes - fesom_salt_first_nodes # Make SideElements - selements_temp = fesom_sidegrid(elements, fesom_temp_last_nodes, lon0, lat_max) - selements_salt = fesom_sidegrid(elements, fesom_salt_last_nodes, lon0, lat_max) - selements_temp_anom = fesom_sidegrid(elements, fesom_temp_anom_nodes, lon0, lat_max) - selements_salt_anom = fesom_sidegrid(elements, fesom_salt_anom_nodes, lon0, lat_max) + selements_temp = fesom_sidegrid(elements, fesom_temp_nodes, lon0, lat_max) + selements_salt = fesom_sidegrid(elements, fesom_salt_nodes, lon0, lat_max) # Build an array of quadrilateral patches for the plot, and of data values # corresponding to each SideElement patches = [] @@ -121,59 +116,26 @@ def mip_drift_slices (roms_grid, roms_file_first, roms_file_last, fesom_mesh_pat fesom_salt = [] for selm in selements_salt: fesom_salt.append(selm.var) - fesom_temp_anom = [] - for selm in selements_temp_anom: - fesom_temp_anom.append(selm.var) - fesom_salt_anom = [] - for selm in selements_salt_anom: - fesom_salt_anom.append(selm.var) # Set up axis labels the way we want them lat_ticks = arange(lat_min+3, lat_max+10, 10) lat_labels = [] for val in lat_ticks: lat_labels.append(str(int(round(-val))) + r'$^{\circ}$S') - depth_ticks = range(depth_min+300, 0+1000, 1000) + depth_ticks = range(depth_min+1000, 0+1000, 1000) depth_labels = [] for val in depth_ticks: depth_labels.append(str(int(round(-val)))) - + print 'Plotting' - fig = figure(figsize=(16,25)) - # Absolute values - gs1 = GridSpec(2,2) - gs1.update(left=0.1, right=0.95, bottom=0.57, top=0.925, wspace=0.08, hspace=0.2) - # MetROMS temperature + fig = figure(figsize=(14,18)) + # ECCO2 + gs1 = GridSpec(1,2) + gs1.update(left=0.1, right=0.95, bottom=0.69, top=0.93, wspace=0.08) + # Temperature ax = subplot(gs1[0,0]) - pcolor(roms_lat, roms_z, roms_temp, vmin=temp_min, vmax=temp_max, cmap='jet') - title(r'MetROMS temperature ($^{\circ}$)', fontsize=24) - ylabel('Depth (m)', fontsize=18) - xlim([lat_min, lat_max]) - ylim([depth_min, depth_max]) - ax.set_xticks(lat_ticks) - ax.set_xticklabels([]) - ax.set_yticks(depth_ticks) - ax.set_yticklabels(depth_labels, fontsize=16) - text(-38, 1000, 'Last year, ' + lon_string, fontsize=30) - # MetROMS salinity - ax = subplot(gs1[0,1]) - pcolor(roms_lat, roms_z, roms_salt, vmin=salt_min, vmax=salt_max, cmap='jet') - title('MetROMS salinity (psu)', fontsize=24) - xlim([lat_min, lat_max]) - ylim([depth_min, depth_max]) - ax.set_xticks(lat_ticks) - ax.set_xticklabels([]) - ax.set_yticks(depth_ticks) - ax.set_yticklabels([]) - # FESOM temperature - ax = subplot(gs1[1,0]) - img = PatchCollection(patches, cmap='jet') - img.set_array(array(fesom_temp)) - 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=24) - xlabel('Latitude', fontsize=18) + pcolor(ecco_lat, ecco_depth, ecco_temp, vmin=temp_min, vmax=temp_max, cmap='jet') + title(r'Temperature ($^{\circ}$C)', fontsize=24) ylabel('Depth (m)', fontsize=18) xlim([lat_min, lat_max]) ylim([depth_min, depth_max]) @@ -181,63 +143,50 @@ def mip_drift_slices (roms_grid, roms_file_first, roms_file_last, fesom_mesh_pat ax.set_xticklabels(lat_labels, fontsize=16) ax.set_yticks(depth_ticks) ax.set_yticklabels(depth_labels, fontsize=16) - # Add a colorbar for temperature - cbaxes = fig.add_axes([0.17, 0.51, 0.3, 0.02]) - cbar = colorbar(img, orientation='horizontal', cax=cbaxes, extend='both', ticks=arange(temp_min, temp_max+2, 2)) - cbar.ax.tick_params(labelsize=16) - # FESOM salinity - ax = subplot(gs1[1,1]) - img = PatchCollection(patches, cmap='jet') - img.set_array(array(fesom_salt)) - img.set_edgecolor('face') - img.set_clim(vmin=salt_min, vmax=salt_max) - ax.add_collection(img) - title('FESOM (high-res) salinity (psu)', fontsize=24) - xlabel('Latitude', fontsize=18) + text(-54, 1000, 'a) ECCO2 initial conditions at ' + lon_string, fontsize=28) + # Salinity + ax = subplot(gs1[0,1]) + pcolor(ecco_lat, ecco_depth, ecco_salt, vmin=salt_min, vmax=salt_max, cmap='jet') + title('Salinity (psu)', fontsize=24) xlim([lat_min, lat_max]) ylim([depth_min, depth_max]) ax.set_xticks(lat_ticks) ax.set_xticklabels(lat_labels, fontsize=16) ax.set_yticks(depth_ticks) ax.set_yticklabels([]) - # Add a colorbar for salinity - cbaxes = fig.add_axes([0.6, 0.51, 0.3, 0.02]) - cbar = colorbar(img, orientation='horizontal', cax=cbaxes, extend='both', ticks=arange(salt_min+0.1, salt_max+0.1, 0.2)) - cbar.ax.tick_params(labelsize=16) - # Anomalies for last year - gs2 = GridSpec(2,2) - gs2.update(left=0.1, right=0.95, bottom=0.075, top=0.43, wspace=0.08, hspace=0.2) - # MetROMS temperature + # MetROMS + gs2 = GridSpec(1,2) + gs2.update(left=0.1, right=0.95, bottom=0.38, top=0.62, wspace=0.08) + # Temperature ax = subplot(gs2[0,0]) - pcolor(roms_lat, roms_z, roms_temp_anom, vmin=-temp_anom, vmax=temp_anom, cmap='RdBu_r') - title(r'MetROMS temperature ($^{\circ}$)', fontsize=24) + pcolor(roms_lat, roms_z, roms_temp, vmin=temp_min, vmax=temp_max, cmap='jet') ylabel('Depth (m)', fontsize=18) xlim([lat_min, lat_max]) ylim([depth_min, depth_max]) ax.set_xticks(lat_ticks) - ax.set_xticklabels([]) + ax.set_xticklabels(lat_labels, fontsize=16) ax.set_yticks(depth_ticks) ax.set_yticklabels(depth_labels, fontsize=16) - text(-47, 1000, 'Last year minus first year, ' + lon_string, fontsize=30) - # MetROMS salinity + text(-55, 300, 'b) MetROMS, last 5 days of simulation', fontsize=28) + # Salinity ax = subplot(gs2[0,1]) - pcolor(roms_lat, roms_z, roms_salt_anom, vmin=-salt_anom, vmax=salt_anom, cmap='RdBu_r') - title('MetROMS salinity (psu)', fontsize=24) + pcolor(roms_lat, roms_z, roms_salt, vmin=salt_min, vmax=salt_max, cmap='jet') xlim([lat_min, lat_max]) ylim([depth_min, depth_max]) ax.set_xticks(lat_ticks) - ax.set_xticklabels([]) + ax.set_xticklabels(lat_labels, fontsize=16) ax.set_yticks(depth_ticks) ax.set_yticklabels([]) - # FESOM temperature - ax = subplot(gs2[1,0]) - img = PatchCollection(patches, cmap='RdBu_r') - img.set_array(array(fesom_temp_anom)) + # FESOM + gs3 = GridSpec(1,2) + gs3.update(left=0.1, right=0.95, bottom=0.07, top=0.31, wspace=0.08) + # Temperature + ax = subplot(gs3[0,0]) + img = PatchCollection(patches, cmap='jet') + img.set_array(array(fesom_temp)) img.set_edgecolor('face') - img.set_clim(vmin=-temp_anom, vmax=temp_anom) + img.set_clim(vmin=temp_min, vmax=temp_max) ax.add_collection(img) - title(r'FESOM (high-res) temperature ($^{\circ}$C)', fontsize=24) - xlabel('Latitude', fontsize=18) ylabel('Depth (m)', fontsize=18) xlim([lat_min, lat_max]) ylim([depth_min, depth_max]) @@ -245,19 +194,18 @@ def mip_drift_slices (roms_grid, roms_file_first, roms_file_last, fesom_mesh_pat ax.set_xticklabels(lat_labels, fontsize=16) ax.set_yticks(depth_ticks) ax.set_yticklabels(depth_labels, fontsize=16) + text(-59, 300, 'c) FESOM (high-res), last 5 days of simulation', fontsize=28) # Add a colorbar for temperature cbaxes = fig.add_axes([0.17, 0.02, 0.3, 0.02]) - cbar = colorbar(img, orientation='horizontal', cax=cbaxes, extend='both', ticks=arange(-temp_anom, temp_anom+1, 1)) + cbar = colorbar(img, orientation='horizontal', cax=cbaxes, extend='both', ticks=arange(temp_min, temp_max+2, 2)) cbar.ax.tick_params(labelsize=16) - # FESOM salinity - ax = subplot(gs2[1,1]) - img = PatchCollection(patches, cmap='RdBu_r') - img.set_array(array(fesom_salt_anom)) + # Salinity + ax = subplot(gs3[0,1]) + img = PatchCollection(patches, cmap='jet') + img.set_array(array(fesom_salt)) img.set_edgecolor('face') - img.set_clim(vmin=-salt_anom, vmax=salt_anom) + img.set_clim(vmin=salt_min, vmax=salt_max) ax.add_collection(img) - title('FESOM (high-res) salinity (psu)', fontsize=24) - xlabel('Latitude', fontsize=18) xlim([lat_min, lat_max]) ylim([depth_min, depth_max]) ax.set_xticks(lat_ticks) @@ -266,24 +214,19 @@ def mip_drift_slices (roms_grid, roms_file_first, roms_file_last, fesom_mesh_pat ax.set_yticklabels([]) # Add a colorbar for salinity cbaxes = fig.add_axes([0.6, 0.02, 0.3, 0.02]) - cbar = colorbar(img, orientation='horizontal', cax=cbaxes, extend='both', ticks=arange(-salt_anom, salt_anom+0.25, 0.25)) + cbar = colorbar(img, orientation='horizontal', cax=cbaxes, extend='both', ticks=arange(salt_min+0.1, salt_max+0.1, 0.2)) cbar.ax.tick_params(labelsize=16) #fig.show() - fig.savefig('ts_slices.png') + fig.savefig('ts_drift.png') # Command-line interface if __name__ == "__main__": roms_grid = raw_input("Path to ROMS grid file: ") - roms_file_first = raw_input("Path to ROMS temperature/salinity file averaged over first year: ") - roms_file_last = raw_input("Path to ROMS temperature/salinity file averaged over last year: ") + roms_file = raw_input("Path to last ocean_avg file output by ROMS: ") fesom_mesh_path = raw_input("Path to FESOM mesh directory: ") - fesom_file_first = raw_input("Path to FESOM temperature/salinity file averaged over first year: ") - fesom_file_last = raw_input("Path to FESOM temperature/salinity file averaged over last year: ") - mip_drift_slices(roms_grid, roms_file_first, roms_file_last, fesom_mesh_path, fesom_file_first, fesom_file_last) - - + fesom_file = raw_input("Path to last oce.mean file output by FESOM: ") + mip_drift_slices(roms_grid, roms_file, fesom_mesh_path, fesom_file) - diff --git a/mip_massloss_difference_map.py b/mip_massloss_difference_map.py index f3ac776..35cef85 100644 --- a/mip_massloss_difference_map.py +++ b/mip_massloss_difference_map.py @@ -5,17 +5,19 @@ from matplotlib.colors import LinearSegmentedColormap # Make a circumpolar Antarctic figure showing the percentage change in mass -# loss for each ice shelf in the FESOM simulation with respect to the MetROMS -# simulation, discarding the first 10 years (i.e. 2002-2016 average). +# loss for each ice shelf in the high-resolution FESOM simulation with respect +# to the low-resolution FESOM simulation, discarding the first 10 years +# (i.e. 2002-2016 average). Also print the values to the screen. # Input: # roms_grid = path to ROMS grid file -# roms_logfile = path to ROMS logfile from timeseries_massloss.py -# fesom_logfile = path to FESOM logfile from timeseries_massloss.py in the -# fesomtools repository +# fesom_logfile_lr, fesom_logfile_hr = path to FESOM logfiles from +# timeseries_massloss.py in the fesomtools repository, for +# the low-resolution and high-resolution simulations +# respectively # save = optional boolean indicating to save the figure to a file, rather than # display on screen # fig_name = if save=True, filename for figure -def mip_massloss_difference_map (roms_grid, roms_logfile, fesom_logfile, save=False, fig_name=None): +def mip_massloss_difference_map (roms_grid, fesom_logfile_lr, fesom_logfile_hr, save=False, fig_name=None): # Year simulations start year_start = 1992 @@ -33,6 +35,8 @@ def mip_massloss_difference_map (roms_grid, roms_logfile, fesom_logfile, save=Fa lon_max = [-59.33, -60, -66.67, -28.33, -88.83, -99.17, -103.33, -111.5, -114.33, -140, -145, 146.62, 123.33, 102.5, 89.17, 75, 37.67, 33.33, 16.17, 12.88, 7.6, -10.33, -146.67, 181] lat_min = [-73.03, -69.35, -74.17, -83.5, -73.28, -75.5, -75.5, -75.33, -74.9, -76.42, -78, -67.83, -67.17, -66.67, -67.83, -73.67, -69.83, -71.67, -70.5, -70.75, -71.83, -76.33, -85, -84.5] lat_max = [-69.37, -66.13, -69.5, -74.67, -71.67, -74.17, -74.67, -73.67, -73, -75.17, -76.41, -66.67, -66.5, -64.83, -66.17, -68.33, -68.67, -68.33, -69.33, -69.83, -69.33, -71.5, -77.77, -77] + # Name of each ice 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'] num_shelves = len(lon_min)-1 # Degrees to radians conversion factor @@ -47,78 +51,61 @@ def mip_massloss_difference_map (roms_grid, roms_logfile, fesom_logfile, save=Fa # Minimum zice in ROMS grid min_zice = -10 - # Read ROMS logfile - roms_time = [] - f = open(roms_logfile, 'r') - # Skip the first line (header for time array) + # Read FESOM timeseries + # Start with low-res + tmp = [] + f = open(fesom_logfile_lr, 'r') + # Skip the first line (header) f.readline() - for line in f: - try: - roms_time.append(float(line)) - except(ValueError): - # Reached the header for the next variable - break - # Skip the values for the entire continent + # Count the number of time indices for the first variable (total mass loss + # for all ice shelves, which we don't care about) + num_time = 0 for line in f: try: tmp = float(line) + num_time += 1 except(ValueError): + # Reached the header for the next variable break # Set up array for mass loss values at each ice shelf - roms_massloss_ts = empty([num_shelves, len(roms_time)]) - index = 0 + fesom_massloss_ts_lr = empty([num_shelves, num_time]) # Loop over ice shelves + index = 0 while index < num_shelves: t = 0 for line in f: try: - roms_massloss_ts[index, t] = float(line) + fesom_massloss_ts_lr[index,t] = float(line) t += 1 except(ValueError): # Reached the header for the next ice shelf break - index +=1 - # Add start year to ROMS time array - roms_time = array(roms_time) + year_start + index += 1 # Average between observation years - t_start = nonzero(roms_time >= obs_start)[0][0] - if obs_end == 2016: - t_end = size(roms_time) - else: - t_end = nonzero(roms_time >= obs_end+1)[0][0] - roms_massloss = mean(roms_massloss_ts[:,t_start:t_end], axis=1) - - # Read FESOM timeseries + fesom_massloss_lr = mean(fesom_massloss_ts_lr[:,peryear*(obs_start-year_start):peryear*(obs_end+1-year_start)], axis=1) + # Repeat for high-res tmp = [] - f = open(fesom_logfile, 'r') - # Skip the first line (header) + f = open(fesom_logfile_hr, 'r') f.readline() - # Count the number of time indices for the first variable (total mass loss - # for all ice shelves, which we don't care about) num_time = 0 for line in f: try: tmp = float(line) num_time += 1 except(ValueError): - # Reached the header for the next variable break - # Set up array for mass loss values at each ice shelf - fesom_massloss_ts = empty([num_shelves, num_time]) - # Loop over ice shelves + fesom_massloss_ts_hr = empty([num_shelves, num_time]) index = 0 while index < num_shelves: t = 0 for line in f: try: - fesom_massloss_ts[index,t] = float(line) + fesom_massloss_ts_hr[index,t] = float(line) t += 1 except(ValueError): - # Reached the header for the next ice shelf break index += 1 - # Average between observation years - fesom_massloss = mean(fesom_massloss_ts[:,peryear*(obs_start-year_start):peryear*(obs_end+1-year_start)], axis=1) + fesom_massloss_hr = mean(fesom_massloss_ts_hr[:,peryear*(obs_start-year_start):peryear*(obs_end+1-year_start)], axis=1) # Read ROMS grid id = Dataset(roms_grid, 'r') @@ -141,8 +128,9 @@ def mip_massloss_difference_map (roms_grid, roms_logfile, fesom_logfile, save=Fa massloss_change[:,:] = ma.masked # Loop over ice shelves for index in range(num_shelves): - # Calculate percentage change in massloss, FESOM wrt MetROMS - tmp = (fesom_massloss[index] - roms_massloss[index])/roms_massloss[index]*100 + # Calculate percentage change in massloss, high-res wrt low-res + tmp = (fesom_massloss_hr[index] - fesom_massloss_lr[index])/fesom_massloss_lr[index]*100 + print names[index] + ': ' + str(tmp) # Modify plotting field for this region if index == num_shelves-1: # Ross region is split into two @@ -164,11 +152,11 @@ def mip_massloss_difference_map (roms_grid, roms_logfile, fesom_logfile, save=Fa land_circle = zeros(shape(x_reg)) land_circle = ma.masked_where(sqrt((x_reg-x_c)**2 + (y_reg-y_c)**2) > radius, land_circle) - # Set up colour scale, -100 to 175, 0 is white - bound = 170.0 - min_colour = (-100.0+bound)/(2.0*bound) - max_colour = 1.0 - lev = linspace(-100, bound, num=50) + # Set up colour scale, -100 to max value, 0 is white + bound = amax(abs(massloss_change)) + lev = linspace(amin(massloss_change), amax(massloss_change), num=50) + min_colour = (amin(massloss_change)+bound)/(2.0*bound) + max_colour = 1 new_cmap = truncate_colormap(get_cmap('RdBu_r'), min_colour, max_colour) # Plot @@ -182,14 +170,14 @@ def mip_massloss_difference_map (roms_grid, roms_logfile, fesom_logfile, save=Fa contourf(x_reg, y_reg, land_circle, 1, colors=(('0.6', '0.6', '0.6'))) # Now shade the percentage change in mass loss contourf(x, y, massloss_change, lev, cmap=new_cmap) - cbar = colorbar(ticks=arange(-100, bound, 50)) + cbar = colorbar(ticks=arange(0, bound+25, 25)) cbar.ax.tick_params(labelsize=20) # Add a black contour for the ice shelf front rcParams['contour.negative_linestyle'] = 'solid' contour(x, y, zice, levels=[min_zice], colors=('black')) xlim([-nbdry, nbdry]) ylim([-nbdry, nbdry]) - title('% Change in Ice Shelf Mass Loss ('+str(obs_start)+'-'+str(obs_end)+' average)\nFESOM with respect to MetROMS', fontsize=30) + title('% Change in Ice Shelf Mass Loss ('+str(obs_start)+'-'+str(obs_end)+' average)\nfrom increased resolution in FESOM', fontsize=30) axis('off') # Finished @@ -210,8 +198,8 @@ def truncate_colormap(cmap, minval=0.0, maxval=1.0, n=-1): if __name__ == "__main__": roms_grid = raw_input("Path to ROMS grid file: ") - roms_logfile = raw_input("Path to ROMS mass loss logfile: ") - fesom_logfile = raw_input("Path to FESOM mass loss logfile: ") + fesom_logfile_lr = raw_input("Path to FESOM low-res mass loss logfile: ") + fesom_logfile_hr = raw_input("Path to FESOM high-res mass loss logfile: ") action = raw_input("Save figure (s) or display in window (d)? ") if action == 's': save = True @@ -219,5 +207,5 @@ def truncate_colormap(cmap, minval=0.0, maxval=1.0, n=-1): elif action == 'd': save = False fig_name = None - mip_massloss_difference_map(roms_grid, roms_logfile, fesom_logfile, save, fig_name) + mip_massloss_difference_map(roms_grid, fesom_logfile_lr, fesom_logfile_hr, save, fig_name) diff --git a/mip_massloss_error_map.py b/mip_massloss_error_map.py index b04455f..2df5c7e 100644 --- a/mip_massloss_error_map.py +++ b/mip_massloss_error_map.py @@ -3,18 +3,20 @@ from matplotlib.pyplot import * from matplotlib import rcParams -# Make a 2x1 circumpolar Antarctic map of percentage error in mass loss for -# each ice shelf, outside the bounds given by Rignot et al. (2013), in MetROMS -# (left) and FESOM (right), for the 2003-2008 average. +# Make a 3x1 circumpolar Antarctic map of percentage error in mass loss for +# each ice shelf, outside the bounds given by Rignot et al. (2013), in MetROMS, +# FESOM low-res, and FESOM high-res, for the 2003-2008 average. Also print the +# values to the screen. # Input: # roms_grid = path to ROMS grid file # roms_logfile = path to ROMS logfile from timeseries_massloss.py -# fesom_logfile = path to FESOM logfile from timeseries_massloss.py in the -# fesomtools repository +# fesom_logfile_lr, fesom_logfile_hr = paths to FESOM logfiles from +# timeseries_massloss.py in the fesomtools repository, for +# low-res and high-res respectively # save = optional boolean indicating to save the figure to a file, rather than # display on screen # fig_name = if save=True, filename for figure -def mip_massloss_error_map (roms_grid, roms_logfile, fesom_logfile, save=False, fig_name=None): +def mip_massloss_error_map (roms_grid, roms_logfile, fesom_logfile_lr, fesom_logfile_hr, save=False, fig_name=None): # Year simulations start year_start = 1992 @@ -32,6 +34,8 @@ def mip_massloss_error_map (roms_grid, roms_logfile, fesom_logfile, save=False, lon_max = [-59.33, -60, -66.67, -28.33, -88.83, -99.17, -103.33, -111.5, -114.33, -140, -145, 146.62, 123.33, 102.5, 89.17, 75, 37.67, 33.33, 16.17, 12.88, 7.6, -10.33, -146.67, 181] lat_min = [-73.03, -69.35, -74.17, -83.5, -73.28, -75.5, -75.5, -75.33, -74.9, -76.42, -78, -67.83, -67.17, -66.67, -67.83, -73.67, -69.83, -71.67, -70.5, -70.75, -71.83, -76.33, -85, -84.5] lat_max = [-69.37, -66.13, -69.5, -74.67, -71.67, -74.17, -74.67, -73.67, -73, -75.17, -76.41, -66.67, -66.5, -64.83, -66.17, -68.33, -68.67, -68.33, -69.33, -69.83, -69.33, -71.5, -77.77, -77] + # Name of each ice 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'] # 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] @@ -88,8 +92,9 @@ def mip_massloss_error_map (roms_grid, roms_logfile, fesom_logfile, save=False, roms_massloss = mean(roms_massloss_ts[:,t_start:t_end], axis=1) # Read FESOM timeseries + # Low-res tmp = [] - f = open(fesom_logfile, 'r') + f = open(fesom_logfile_lr, 'r') # Skip the first line (header) f.readline() # Count the number of time indices for the first variable (total mass loss @@ -103,21 +108,44 @@ def mip_massloss_error_map (roms_grid, roms_logfile, fesom_logfile, save=False, # Reached the header for the next variable break # Set up array for mass loss values at each ice shelf - fesom_massloss_ts = empty([num_shelves, num_time]) + fesom_massloss_ts_lr = empty([num_shelves, num_time]) # Loop over ice shelves index = 0 while index < num_shelves: t = 0 for line in f: try: - fesom_massloss_ts[index,t] = float(line) + fesom_massloss_ts_lr[index,t] = float(line) t += 1 except(ValueError): # Reached the header for the next ice shelf break index += 1 # Average between observation years - fesom_massloss = mean(fesom_massloss_ts[:,peryear*(obs_start-year_start):peryear*(obs_end+1-year_start)], axis=1) + fesom_massloss_lr = mean(fesom_massloss_ts_lr[:,peryear*(obs_start-year_start):peryear*(obs_end+1-year_start)], axis=1) + # Repeat for high-res + tmp = [] + f = open(fesom_logfile_hr, 'r') + f.readline() + num_time = 0 + for line in f: + try: + tmp = float(line) + num_time += 1 + except(ValueError): + break + fesom_massloss_ts_hr = empty([num_shelves, num_time]) + index = 0 + while index < num_shelves: + t = 0 + for line in f: + try: + fesom_massloss_ts_hr[index,t] = float(line) + t += 1 + except(ValueError): + break + index += 1 + fesom_massloss_hr = mean(fesom_massloss_ts_hr[:,peryear*(obs_start-year_start):peryear*(obs_end+1-year_start)], axis=1) # Read ROMS grid id = Dataset(roms_grid, 'r') @@ -137,11 +165,14 @@ def mip_massloss_error_map (roms_grid, roms_logfile, fesom_logfile, save=False, # Initialise fields of ice shelf mass loss unexplained percent error in each mode roms_error = ma.empty(shape(lon)) - fesom_error = ma.empty(shape(lon)) + fesom_error_lr = ma.empty(shape(lon)) + fesom_error_hr = ma.empty(shape(lon)) roms_error[:,:] = ma.masked - fesom_error[:,:] = ma.masked + fesom_error_lr[:,:] = ma.masked + fesom_error_hr[:,:] = ma.masked # Loop over ice shelves for index in range(num_shelves): + print names[index] # Find the range of observations massloss_low = obs_massloss[index] - obs_massloss_error[index] massloss_high = obs_massloss[index] + obs_massloss_error[index] @@ -156,6 +187,7 @@ def mip_massloss_error_map (roms_grid, roms_logfile, fesom_logfile, save=False, else: # Simulated mass loss within observational error estimates error_tmp = 0 + print 'MetROMS: ' + str(error_tmp) # Modify error field for this region if index == num_shelves-1: # Ross region is split into two @@ -163,23 +195,42 @@ def mip_massloss_error_map (roms_grid, roms_logfile, fesom_logfile, save=False, else: region = (lon >= lon_min[index])*(lon <= lon_max[index])*(lat >= lat_min[index])*(lat <= lat_max[index])*(mask_zice == 1) roms_error[region] = error_tmp - # Repeat for FESOM - if fesom_massloss[index] < massloss_low: + # Repeat for FESOM low-res + if fesom_massloss_lr[index] < massloss_low: # Simulated mass loss too low - error_tmp = (fesom_massloss[index] - massloss_low)/massloss_low*100 - elif fesom_massloss[index] > massloss_high: + error_tmp = (fesom_massloss_lr[index] - massloss_low)/massloss_low*100 + elif fesom_massloss_lr[index] > massloss_high: # Simulated mass loss too high - error_tmp = (fesom_massloss[index] - massloss_high)/massloss_high*100 + error_tmp = (fesom_massloss_lr[index] - massloss_high)/massloss_high*100 else: # Simulated mass loss within observational error estimates error_tmp = 0 + print 'FESOM low-res: ' + str(error_tmp) # Modify error field for this region if index == num_shelves-1: # Ross region is split into two region = (lon >= lon_min[index])*(lon <= lon_max[index])*(lat >= lat_min[index])*(lat <= lat_max[index])*(mask_zice == 1) + (lon >= lon_min[index+1])*(lon <= lon_max[index+1])*(lat >= lat_min[index+1])*(lat <= lat_max[index+1])*(mask_zice == 1) else: region = (lon >= lon_min[index])*(lon <= lon_max[index])*(lat >= lat_min[index])*(lat <= lat_max[index])*(mask_zice == 1) - fesom_error[region] = error_tmp + fesom_error_lr[region] = error_tmp + # Repeat for FESOM high-res + if fesom_massloss_hr[index] < massloss_low: + # Simulated mass loss too low + error_tmp = (fesom_massloss_hr[index] - massloss_low)/massloss_low*100 + elif fesom_massloss_hr[index] > massloss_high: + # Simulated mass loss too high + error_tmp = (fesom_massloss_hr[index] - massloss_high)/massloss_high*100 + else: + # Simulated mass loss within observational error estimates + error_tmp = 0 + print 'FESOM high-res: ' + str(error_tmp) + # Modify error field for this region + if index == num_shelves-1: + # Ross region is split into two + region = (lon >= lon_min[index])*(lon <= lon_max[index])*(lat >= lat_min[index])*(lat <= lat_max[index])*(mask_zice == 1) + (lon >= lon_min[index+1])*(lon <= lon_max[index+1])*(lat >= lat_min[index+1])*(lat <= lat_max[index+1])*(mask_zice == 1) + else: + region = (lon >= lon_min[index])*(lon <= lon_max[index])*(lat >= lat_min[index])*(lat <= lat_max[index])*(mask_zice == 1) + fesom_error_hr[region] = error_tmp # Edit zice so tiny ice shelves won't be contoured zice[roms_error.mask] = 0.0 @@ -197,10 +248,12 @@ def mip_massloss_error_map (roms_grid, roms_logfile, fesom_logfile, save=False, lev = linspace(-100, 100, num=50) # Plot - fig = figure(figsize=(20,9)) + fig = figure(figsize=(19,8)) fig.patch.set_facecolor('white') + gs = GridSpec(1,3) + gs.update(left=0, right=1, bottom=0.1, top=0.85, wspace=0.05) # ROMS - ax1 = fig.add_subplot(1,2,1, aspect='equal') + ax1 = subplot(gs[0,0], aspect='equal') # First shade land and zice in grey (include zice so there are no white # patches near the grounding line where contours meet) contourf(x, y, land_zice, 1, colors=(('0.6', '0.6', '0.6'))) @@ -213,26 +266,36 @@ def mip_massloss_error_map (roms_grid, roms_logfile, fesom_logfile, save=False, contour(x, y, zice, levels=[min_zice], colors=('black')) xlim([-nbdry, nbdry]) ylim([-nbdry, nbdry]) - title('MetROMS', fontsize=30) + title('a) MetROMS', fontsize=28) + axis('off') + # FESOM low-res + ax2 = subplot(gs[0,1], aspect='equal') + contourf(x, y, land_zice, 1, colors=(('0.6', '0.6', '0.6'))) + contourf(x_reg, y_reg, land_circle, 1, colors=(('0.6', '0.6', '0.6'))) + img = contourf(x, y, fesom_error_lr, lev, cmap='RdBu_r', extend='max') + rcParams['contour.negative_linestyle'] = 'solid' + contour(x, y, zice, levels=[min_zice], colors=('black')) + xlim([-nbdry, nbdry]) + ylim([-nbdry, nbdry]) axis('off') - # FESOM - ax2 = fig.add_subplot(1,2,2, aspect='equal') + title('b) FESOM low-res', fontsize=28) + # FESOM high-res + ax3 = subplot(gs[0,2], aspect='equal') contourf(x, y, land_zice, 1, colors=(('0.6', '0.6', '0.6'))) contourf(x_reg, y_reg, land_circle, 1, colors=(('0.6', '0.6', '0.6'))) - img = contourf(x, y, fesom_error, lev, cmap='RdBu_r', extend='max') + img = contourf(x, y, fesom_error_hr, lev, cmap='RdBu_r', extend='max') rcParams['contour.negative_linestyle'] = 'solid' contour(x, y, zice, levels=[min_zice], colors=('black')) xlim([-nbdry, nbdry]) ylim([-nbdry, nbdry]) axis('off') - title('FESOM', fontsize=30) - # Add a colourbar on the right - cbaxes = fig.add_axes([0.92, 0.2, 0.01, 0.6]) - cbar = colorbar(img, cax=cbaxes, ticks=arange(-100, 100+25, 25)) + title('c) FESOM high-res', fontsize=28) + # Add a horizontal colourbar on the bottom + cbaxes = fig.add_axes([0.34, 0.1, 0.4, 0.04]) + cbar = colorbar(img, orientation='horizontal', cax=cbaxes, ticks=arange(-100, 100+25, 25)) cbar.ax.tick_params(labelsize=20) # Main title - suptitle('Bias in Ice Shelf Mass Loss (%)\n2003-2008', fontsize=24) - subplots_adjust(wspace=0.05) + suptitle('Bias in Ice Shelf Mass Loss (%), 2003-2008', fontsize=34) # Finished if save: @@ -241,11 +304,13 @@ def mip_massloss_error_map (roms_grid, roms_logfile, fesom_logfile, save=False, fig.show() +# Command-line interface if __name__ == "__main__": roms_grid = raw_input("Path to ROMS grid file: ") roms_logfile = raw_input("Path to ROMS mass loss logfile: ") - fesom_logfile = raw_input("Path to FESOM mass loss logfile: ") + fesom_logfile_lr = raw_input("Path to FESOM low-res mass loss logfile: ") + fesom_logfile_hr = raw_input("Path to FESOM high-res mass loss logfile: ") action = raw_input("Save figure (s) or display in window (d)? ") if action == 's': save = True @@ -253,5 +318,5 @@ def mip_massloss_error_map (roms_grid, roms_logfile, fesom_logfile, save=False, elif action == 'd': save = False fig_name = None - mip_massloss_error_map (roms_grid, roms_logfile, fesom_logfile, save, fig_name) + mip_massloss_error_map (roms_grid, roms_logfile, fesom_logfile_lr, fesom_logfile_hr, save, fig_name) diff --git a/mip_ts_distribution.py b/mip_ts_distribution.py index e124d9d..72ae325 100644 --- a/mip_ts_distribution.py +++ b/mip_ts_distribution.py @@ -208,13 +208,13 @@ def mip_ts_distribution (roms_grid, roms_file, fesom_mesh_path, fesom_file): title('FESOM (high-res)', fontsize=28) # Add a colourbar on the right cbaxes = fig.add_axes([0.93, 0.2, 0.02, 0.6]) - cbar = colorbar(img, cax=cbaxes, ticks=[0,100,500,1000,2000,3000,4000]) + cbar = colorbar(img, cax=cbaxes, ticks=[0,50,100,200,500,1000,2000,4000]) cbar.ax.tick_params(labelsize=18) # Add the main title suptitle('Water masses south of 65$^{\circ}$S: depth (m), 2002-2016 average', fontsize=30) subplots_adjust(wspace=0.1) #fig.show() - fig.savefig('ts_distribution.png') + fig.savefig('ts_distribution_orig.png') # Command-line interface diff --git a/zonal_plot.py b/zonal_plot.py index 932c663..344e694 100644 --- a/zonal_plot.py +++ b/zonal_plot.py @@ -79,6 +79,10 @@ def zonal_plot (file_path, var_name, tstep, lon_key, lon0, lon_bounds, depth_min # Get a 3D array of z-coordinates; sc_r and Cs_r are unused in this script z_3d, sc_r, Cs_r = calc_z(h, zice, theta_s, theta_b, hc, N, zeta) + # Warning message for zonal averages + if lon_key != 0: + print 'WARNING: this script assumes regular i-indices for zonal averages. In heavily rotated regions eg inner Weddell Sea it may not be appropriate.' + # Choose what to write on the title about longitude if lon_key == 0: if lon0 < 0: