From 747398a5dc9fb47bf562aca517f98af01c8b86ea Mon Sep 17 00:00:00 2001 From: Kaitlin Naughten Date: Wed, 14 Jun 2017 15:55:59 +1000 Subject: [PATCH] Added intercomparison figures for timeseries and mass loss maps --- file_guide.txt | 55 +++++ massloss_map.py | 2 +- mip_massloss_difference_map.py | 206 +++++++++++++++++++ mip_massloss_error_map.py | 257 ++++++++++++++++++++++++ mip_timeseries.py | 353 +++++++++++++++++++++++++++++++++ timeseries_massloss.py | 1 + 6 files changed, 873 insertions(+), 1 deletion(-) create mode 100644 mip_massloss_difference_map.py create mode 100644 mip_massloss_error_map.py create mode 100644 mip_timeseries.py diff --git a/file_guide.txt b/file_guide.txt index 4b6abbf..9c5eb4e 100644 --- a/file_guide.txt +++ b/file_guide.txt @@ -1115,6 +1115,61 @@ mip_grid_res.py: Make a 2x1 plot showing horizontal grid resolution (square directory, and whether you want to save the figure (and if so, what filename) or display it on the screen. +mip_timeseries.py: Plot MetROMS and FESOM timeseries together, for Drake + Passage transport, total Antarctic sea ice area and volume, + and basal mass loss for major ice shelves. Include the range + of observations for Drake Passage transport and ice shelf + mass loss. + To run: First make sure you have run timeseries_dpt.py, + timeseries_seaice.py, and timeseries_massloss.py + and saved the logfiles in a single directory with + the names dpt.log, seaice.log, and massloss.log. Do + the same for FESOM with the equivalent fesomtools + scripts, with the same names. Then open python or + ipython and type "run mip_timeseries.py". The script + will prompt you for the paths to the ROMS and FESOM + logfile directories, and whether you want to plot + the full timeseries or annual averages (in which case + sea ice area and volume will not be plotted). The + script will create a bunch of png files. + +mip_massloss_error_map.py: 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. + 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 + "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. + +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. + 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 + "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. + ***UTILITY FUNCTIONS*** diff --git a/massloss_map.py b/massloss_map.py index 492fc99..22e3ae9 100644 --- a/massloss_map.py +++ b/massloss_map.py @@ -124,7 +124,7 @@ def massloss_map (grid_path, log_path, save=False, fig_name=None): error[region] = error_tmp # Edit zice so tiny ice shelves won't be contoured - #zice[error.mask] = 0.0 + zice[error.mask] = 0.0 # Convert grid to spherical coordinates x = -(lat+90)*cos(lon*deg2rad+pi/2) diff --git a/mip_massloss_difference_map.py b/mip_massloss_difference_map.py new file mode 100644 index 0000000..287538b --- /dev/null +++ b/mip_massloss_difference_map.py @@ -0,0 +1,206 @@ +from netCDF4 import Dataset +from numpy import * +from matplotlib.pyplot import * +from matplotlib import rcParams + +# 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. +# 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 +# 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): + + # Year simulations start + year_start = 1992 + # Years to analyse + obs_start = 2003 + obs_end = 2008 + # Number of output steps per year in FESOM + peryear = 365/5 + + # Limits on longitude and latitude for each ice shelf + # These depend on the source geometry, in this case RTopo 1.05 + # Note there is one extra index at the end of each array; this is because + # the Ross region crosses the line 180W and therefore is split into two + lon_min = [-62.67, -65.5, -79.17, -85, -104.17, -102.5, -108.33, -114.5, -135.67, -149.17, -155, 144, 115, 94.17, 80.83, 65, 33.83, 19, 12.9, 9.33, -10.05, -28.33, -181, 158.33] + 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] + num_shelves = len(lon_min)-1 + + # Degrees to radians conversion factor + deg2rad = pi/180 + # Northern boundary 63S for plot + nbdry = -63+90 + # Centre of missing circle in grid + lon_c = 50 + lat_c = -83 + # Radius of missing circle (play around with this until it works) + radius = 10.1 + # 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) + 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 + for line in f: + try: + tmp = float(line) + except(ValueError): + 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 + # Add start year to ROMS time array + roms_time = array(roms_time) + year_start + # Average between observation years + t_start = nonzero(roms_time >= obs_start)[0][0] + 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 + tmp = [] + f = open(fesom_logfile, 'r') + # Skip the first line (header) + 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 + index = 0 + while index < num_shelves: + t = 0 + for line in f: + try: + fesom_massloss_ts[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) + + # Read ROMS grid + id = Dataset(roms_grid, 'r') + lon = id.variables['lon_rho'][:-15,:-1] + lat = id.variables['lat_rho'][:-15,:-1] + mask_rho = id.variables['mask_rho'][:-15,:-1] + mask_zice = id.variables['mask_zice'][:-15,:-1] + zice = id.variables['zice'][:-15,:-1] + id.close() + # Make sure longitude goes from -180 to 180, not 0 to 360 + index = lon > 180 + lon[index] = lon[index] - 360 + # Get land/zice mask + open_ocn = copy(mask_rho) + open_ocn[mask_zice==1] = 0 + land_zice = ma.masked_where(open_ocn==1, open_ocn) + + # Initialise plotting field of ice shelf mass loss percentage change + massloss_change = ma.empty(shape(lon)) + 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 + # Modify plotting 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) + massloss_change[region] = tmp + + # Edit zice so tiny ice shelves won't be contoured + zice[massloss_change.mask] = 0.0 + # Convert grid to spherical coordinates + x = -(lat+90)*cos(lon*deg2rad+pi/2) + y = (lat+90)*sin(lon*deg2rad+pi/2) + # Find centre in spherical coordinates + x_c = -(lat_c+90)*cos(lon_c*deg2rad+pi/2) + y_c = (lat_c+90)*sin(lon_c*deg2rad+pi/2) + # Build a regular x-y grid and select the missing circle + x_reg, y_reg = meshgrid(linspace(-nbdry, nbdry, num=1000), linspace(-nbdry, nbdry, num=1000)) + 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 + lev = linspace(-200, 200, num=50) + + # Plot + fig = figure(figsize=(16,12)) + ax = fig.add_subplot(1,1,1, aspect='equal') + fig.patch.set_facecolor('white') + # 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'))) + # Fill in the missing cicle + 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='RdBu_r', extend='both') + cbar = colorbar(ticks=arange(-200, 200+50, 50)) + 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 (2003-2008 average)\nFESOM with respect to MetROMS', fontsize=30) + axis('off') + + # Finished + if save: + fig.savefig(fig_name) + else: + fig.show() + + +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: ") + action = raw_input("Save figure (s) or display in window (d)? ") + if action == 's': + save = True + fig_name = raw_input("File name for figure: ") + elif action == 'd': + save = False + fig_name = None + mip_massloss_difference (roms_grid, roms_logfile, fesom_logfile, save, fig_name) + diff --git a/mip_massloss_error_map.py b/mip_massloss_error_map.py new file mode 100644 index 0000000..b04455f --- /dev/null +++ b/mip_massloss_error_map.py @@ -0,0 +1,257 @@ +from netCDF4 import Dataset +from numpy import * +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. +# 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 +# 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): + + # Year simulations start + year_start = 1992 + # Years observations span + obs_start = 2003 + obs_end = 2008 + # Number of output steps per year in FESOM + peryear = 365/5 + + # Limits on longitude and latitude for each ice shelf + # These depend on the source geometry, in this case RTopo 1.05 + # Note there is one extra index at the end of each array; this is because + # the Ross region crosses the line 180W and therefore is split into two + lon_min = [-62.67, -65.5, -79.17, -85, -104.17, -102.5, -108.33, -114.5, -135.67, -149.17, -155, 144, 115, 94.17, 80.83, 65, 33.83, 19, 12.9, 9.33, -10.05, -28.33, -181, 158.33] + 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] + # 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] + num_shelves = len(obs_massloss) + + # Degrees to radians conversion factor + deg2rad = pi/180 + # Northern boundary 63S for plot + nbdry = -63+90 + # Centre of missing circle in grid + lon_c = 50 + lat_c = -83 + # Radius of missing circle (play around with this until it works) + radius = 10.1 + # 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) + 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 + for line in f: + try: + tmp = float(line) + except(ValueError): + 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 + # Add start year to ROMS time array + roms_time = array(roms_time) + year_start + # Average between observation years + t_start = nonzero(roms_time >= obs_start)[0][0] + 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 + tmp = [] + f = open(fesom_logfile, 'r') + # Skip the first line (header) + 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 + index = 0 + while index < num_shelves: + t = 0 + for line in f: + try: + fesom_massloss_ts[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) + + # Read ROMS grid + id = Dataset(roms_grid, 'r') + lon = id.variables['lon_rho'][:-15,:-1] + lat = id.variables['lat_rho'][:-15,:-1] + mask_rho = id.variables['mask_rho'][:-15,:-1] + mask_zice = id.variables['mask_zice'][:-15,:-1] + zice = id.variables['zice'][:-15,:-1] + id.close() + # Make sure longitude goes from -180 to 180, not 0 to 360 + index = lon > 180 + lon[index] = lon[index] - 360 + # Get land/zice mask + open_ocn = copy(mask_rho) + open_ocn[mask_zice==1] = 0 + land_zice = ma.masked_where(open_ocn==1, open_ocn) + + # 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)) + roms_error[:,:] = ma.masked + fesom_error[:,:] = ma.masked + # Loop over ice shelves + for index in range(num_shelves): + # Find the range of observations + massloss_low = obs_massloss[index] - obs_massloss_error[index] + massloss_high = obs_massloss[index] + obs_massloss_error[index] + # Find the unexplained percent error in mass loss + # ROMS + if roms_massloss[index] < massloss_low: + # Simulated mass loss too low + error_tmp = (roms_massloss[index] - massloss_low)/massloss_low*100 + elif roms_massloss[index] > massloss_high: + # Simulated mass loss too high + error_tmp = (roms_massloss[index] - massloss_high)/massloss_high*100 + else: + # Simulated mass loss within observational error estimates + error_tmp = 0 + # 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) + roms_error[region] = error_tmp + # Repeat for FESOM + if fesom_massloss[index] < massloss_low: + # Simulated mass loss too low + error_tmp = (fesom_massloss[index] - massloss_low)/massloss_low*100 + elif fesom_massloss[index] > massloss_high: + # Simulated mass loss too high + error_tmp = (fesom_massloss[index] - massloss_high)/massloss_high*100 + else: + # Simulated mass loss within observational error estimates + error_tmp = 0 + # 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 + + # Edit zice so tiny ice shelves won't be contoured + zice[roms_error.mask] = 0.0 + # Convert grid to spherical coordinates + x = -(lat+90)*cos(lon*deg2rad+pi/2) + y = (lat+90)*sin(lon*deg2rad+pi/2) + # Find centre in spherical coordinates + x_c = -(lat_c+90)*cos(lon_c*deg2rad+pi/2) + y_c = (lat_c+90)*sin(lon_c*deg2rad+pi/2) + # Build a regular x-y grid and select the missing circle + x_reg, y_reg = meshgrid(linspace(-nbdry, nbdry, num=1000), linspace(-nbdry, nbdry, num=1000)) + 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 + lev = linspace(-100, 100, num=50) + + # Plot + fig = figure(figsize=(20,9)) + fig.patch.set_facecolor('white') + # ROMS + ax1 = fig.add_subplot(1,2,1, 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'))) + # Fill in the missing cicle + contourf(x_reg, y_reg, land_circle, 1, colors=(('0.6', '0.6', '0.6'))) + # Now shade the error in mass loss + contourf(x, y, roms_error, lev, cmap='RdBu_r', extend='both') + # 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('MetROMS', fontsize=30) + axis('off') + # FESOM + ax2 = fig.add_subplot(1,2,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') + 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)) + cbar.ax.tick_params(labelsize=20) + # Main title + suptitle('Bias in Ice Shelf Mass Loss (%)\n2003-2008', fontsize=24) + subplots_adjust(wspace=0.05) + + # Finished + if save: + fig.savefig(fig_name) + else: + fig.show() + + +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: ") + action = raw_input("Save figure (s) or display in window (d)? ") + if action == 's': + save = True + fig_name = raw_input("File name for figure: ") + elif action == 'd': + save = False + fig_name = None + mip_massloss_error_map (roms_grid, roms_logfile, fesom_logfile, save, fig_name) + diff --git a/mip_timeseries.py b/mip_timeseries.py new file mode 100644 index 0000000..3b0c24a --- /dev/null +++ b/mip_timeseries.py @@ -0,0 +1,353 @@ +from numpy import * +from matplotlib.pyplot import * + +# Plot MetROMS and FESOM timeseries together, for Drake Passage transport, total +# Antarctic sea ice area and volume, and basal mass loss for major ice shelves. +# Include the range of observations for Drake Passage transport and ice shelf +# mass loss. +# Input: +# roms_dir = path to ROMS directory containing logfiles from timeseries_dpt.py, +# timeseries_seaice.py, and timeseries_massloss.py. It is assumed +# they are saved with the filenames dpt.log, seaice.log, and +# massloss.log. +# fesom_dir = path to FESOM directory containing logfiles from the equivalent +# "fesomtools" scripts with the same names. It is assumed they are +# also saved with the filenames dpt.log, seaice.log, and +# massloss.log. +# annual = optional boolean indicating to calculate annual averages of Drake +# Passage transport and ice shelf mass loss, and to not bother with +# sea ice area and volume. +def mip_timeseries (roms_dir, fesom_dir, annual=False): + + # Averaging period (days) + days_per_output = 5 + year_start = 1992 + year_end = 2016 + # Titles for plotting + model_titles = ['MetROMS', 'FESOM'] + # Colours for plotting + model_colours = ['blue', 'green'] + + # Titles for each ice shelf + shelf_names = ['All Ice Shelves', 'Larsen D Ice Shelf', 'Larsen C Ice Shelf', 'Wilkins & George VI & Stange Ice Shelves', 'Ronne-Filchner Ice Shelf', 'Abbot Ice Shelf', 'Pine Island Glacier Ice Shelf', 'Thwaites Ice Shelf', 'Dotson Ice Shelf', 'Getz Ice Shelf', 'Nickerson Ice Shelf', 'Sulzberger Ice Shelf', 'Mertz Ice Shelf', 'Totten & Moscow University Ice Shelves', 'Shackleton Ice Shelf', 'West Ice Shelf', 'Amery Ice Shelf', 'Prince Harald Ice Shelf', 'Baudouin & Borchgrevink Ice Shelves', 'Lazarev Ice Shelf', 'Nivl Ice Shelf', 'Fimbul & Jelbart & Ekstrom Ice Shelves', 'Brunt & Riiser-Larsen Ice Shelves', 'Ross Ice Shelf'] + # Beginnings of figure names for each ice shelf + fig_names = ['total_massloss', 'larsen_d', 'larsen_c', 'wilkins_georgevi_stange', 'ronne_filchner', 'abbot', 'pig', 'thwaites', 'dotson', 'getz', 'nickerson', 'sulzberger', 'mertz', 'totten_moscowuni', 'shackleton', 'west', 'amery', 'princeharald', 'baudouin_borchgrevink', 'lazarev', 'nivl', 'fimbul_jelbart_ekstrom', 'brunt_riiserlarsen', 'ross'] + num_shelves = len(shelf_names) + + # Bounds of observations for Drake Passage transport + dpt_low = 107 # Lower bound of range of estimates from Cunningham et al 2003, doi:10.1029/2001JC001147: 134 +/- 15-27 Sv + dpt_high = 173.3 # Donohue et al 2016, doi:10.1002/2016GL070319 + # 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] + + # Bounds for plotting Drake Passage transport + dpt_bounds = [80, 200] + # Bounds for plotting each ice shelf mass loss + massloss_bounds_low = [0, -15, -100, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -6, 0, 0, 0, 0, -20, 0] + massloss_bounds_high = [2000, 30, 300, 250, 600, 100, 120, 110, 55, 250, 35, 50, 30, 100, 120, 75, 850, 25, 150, 40, 25, 120, 175, 500] + + # Drake Passage transport + # Read ROMS timeseries + roms_time = [] + roms_dpt = [] + f = open(roms_dir + 'dpt.log', 'r') + # Skip 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 + for line in f: + roms_dpt.append(float(line)) + f.close() + # Add start year to ROMS time array + roms_time = array(roms_time) + year_start + # Read FESOM timeseries + fesom_dpt = [] + f = open(fesom_dir + 'dpt.log', 'r') + # Skip header + f.readline() + for line in f: + fesom_dpt.append(float(line)) + f.close() + # Make FESOM time array (note that FESOM neglects leap years in its output) + fesom_time = arange(len(fesom_dpt))*days_per_output/365. + year_start + if annual: + # Annually average ROMS data + roms_dpt_avg = roms_annual_avg(roms_time, roms_dpt, year_start) + # Annually average FESOM data + fesom_dpt_avg = fesom_annual_avg(fesom_dpt, days_per_output) + # Make new time arrays + roms_time_annual = arange(len(roms_dpt_avg)) + year_start + fesom_time_annual = arange(len(fesom_dpt_avg)) + year_start + # Find last timestep + max_time = max(roms_time_annual[-1], fesom_time_annual[-1]) + else: + max_time = max(roms_time[-1], fesom_time[-1]) + # Plot + fig, ax = subplots(figsize=(10,6)) + if annual: + ax.plot(roms_time_annual, roms_dpt_avg, label=model_titles[0], color=model_colours[0], linewidth=2) + ax.plot(fesom_time_annual, fesom_dpt_avg, label=model_titles[1], color=model_colours[1], linewidth=2) + else: + ax.plot(roms_time, roms_dpt, label=model_titles[0], color=model_colours[0], linewidth=1) + ax.plot(fesom_time, fesom_dpt, label=model_titles[1], color=model_colours[1], linewidth=1) + # Add lines for range of observations + ax.axhline(dpt_low, color='red', linestyle='dashed', linewidth=2, label='observations') + ax.axhline(dpt_high, color='red', linewidth=2, linestyle='dashed') + if annual: + title('Drake Passage Transport (annually averaged)', fontsize=18) + else: + title('Drake Passage Transport', fontsize=18) + xlabel('Year', fontsize=14) + ylabel('Sv', fontsize=14) + xlim([year_start, max_time]) + ylim([dpt_bounds[0], dpt_bounds[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)) + if annual: + fig.savefig('drakepsgtrans_avg.png') + else: + fig.savefig('drakepsgtrans.png') + + # Sea ice area and volume + if not annual: + # Read ROMS timeseries + roms_icearea = [] + roms_icevolume = [] + f = open(roms_dir + 'seaice.log', 'r') + # Skip the first line (header for time array) + f.readline() + for line in f: + # Skip the time values (already have them from DPT) + try: + tmp = float(line) + except(ValueError): + # Reached the header for the next variable + break + # Sea ice area + for line in f: + try: + roms_icearea.append(float(line)) + except(ValueError): + break + # Sea ice volume + roms_icevolume.append(float(line)) + f.close() + # Read FESOM timeseries + fesom_icearea = [] + fesom_icevolume = [] + f = open(fesom_dir + 'seaice.log', 'r') + f.readline() + for line in f: + try: + fesom_icearea.append(float(line)) + except(ValueError): + break + for line in f: + fesom_icevolume.append(float(line)) + f.close() + # Plot sea ice area + fig, ax = subplots(figsize=(10,6)) + ax.plot(roms_time, roms_icearea, label=model_titles[0], color=model_colours[0], linewidth=1) + ax.plot(fesom_time, fesom_icearea, label=model_titles[1], color=model_colours[1], linewidth=1) + title('Antarctic Sea Ice Area', fontsize=18) + xlabel('Year', fontsize=14) + ylabel(r'million km$^2$', fontsize=14) + xlim([year_start, max_time]) + grid(True) + box = ax.get_position() + ax.set_position([box.x0, box.y0, box.width*0.8, box.height]) + ax.legend(loc='center left', bbox_to_anchor=(1,0.5)) + fig.savefig('seaice_area.png') + # Plot sea ice volume + fig, ax = subplots(figsize=(10,6)) + ax.plot(roms_time, roms_icevolume, label=model_titles[0], color=model_colours[0], linewidth=1) + ax.plot(fesom_time, fesom_icevolume, label=model_titles[1], color=model_colours[1], linewidth=1) + title('Antarctic Sea Ice Volume', fontsize=18) + xlabel('Year', fontsize=14) + ylabel(r'million km$^3$', fontsize=14) + xlim([year_start, max_time]) + grid(True) + box = ax.get_position() + ax.set_position([box.x0, box.y0, box.width*0.8, box.height]) + ax.legend(loc='center left', bbox_to_anchor=(1,0.5)) + fig.savefig('seaice_volume.png') + + # Ice shelf mass loss + # Read ROMS timeseries + roms_massloss = empty([num_shelves, size(roms_time)]) + f = open(roms_dir + 'massloss.log', 'r') + # Skip the first line (header for time array) + f.readline() + for line in f: + # Skip the time values (already have them from DPT) + try: + tmp = float(line) + except(ValueError): + # Reached the header for the next variable + break + # Loop over ice shelves + index = 0 + while index < num_shelves: + t = 0 + for line in f: + try: + roms_massloss[index, t] = float(line) + t += 1 + except(ValueError): + # Reached the header for the next ice shelf + break + index += 1 + f.close() + # Read FESOM timeseries + fesom_massloss = empty([num_shelves, size(fesom_time)]) + f = open(fesom_dir + 'massloss.log', 'r') + f.readline() + # Loop over ice shelves + index = 0 + while index < num_shelves: + t = 0 + for line in f: + try: + fesom_massloss[index, t] = float(line) + t += 1 + except(ValueError): + # Reached the header for the next ice shelf + break + index += 1 + f.close() + if annual: + # Annually average ROMS data for each ice shelf + roms_massloss_avg = empty([num_shelves, size(roms_time_annual)]) + for index in range(num_shelves): + roms_massloss_avg[index,:] = roms_annual_avg(roms_time, roms_massloss[index,:], year_start) + # Annually average FESOM data for each ice shelf + fesom_massloss_avg = empty([num_shelves, size(fesom_time_annual)]) + for index in range(num_shelves): + fesom_massloss_avg[index,:] = fesom_annual_avg(fesom_massloss[index,:], days_per_output) + # One plot for each ice shelf + for index in range(num_shelves): + # Calculate range of observations + massloss_low = obs_massloss[index] - obs_massloss_error[index] + massloss_high = obs_massloss[index] + obs_massloss_error[index] + fig, ax = subplots(figsize=(10,6)) + if annual: + ax.plot(roms_time_annual, roms_massloss_avg[index,:], label=model_titles[0], color=model_colours[0], linewidth=2) + ax.plot(fesom_time_annual, fesom_massloss_avg[index,:], label=model_titles[1], color=model_colours[1], linewidth=2) + else: + ax.plot(roms_time, roms_massloss[index,:], label=model_titles[0], color=model_colours[0], linewidth=1) + ax.plot(fesom_time, fesom_massloss[index,:], label=model_titles[1], color=model_colours[1], linewidth=1) + # Add lines for range of observations + ax.axhline(massloss_low, color='red', linestyle='dashed', linewidth=2, label='observations') + ax.axhline(massloss_high, color='red', linewidth=2, linestyle='dashed') + if annual: + title(shelf_names[index] + '\nBasal Mass Loss (Annually Averaged)', fontsize=18) + else: + title(shelf_names[index] + '\nBasal Mass Loss', fontsize=18) + xlabel('Year', fontsize=14) + ylabel('Gt/y', fontsize=14) + xlim([year_start, max_time]) + if not annual: + ylim([massloss_bounds_low[index], massloss_bounds_high[index]]) + grid(True) + box = ax.get_position() + ax.set_position([box.x0, box.y0, box.width*0.8, box.height]) + ax.legend(loc='center left', bbox_to_anchor=(1,0.5)) + if annual: + fig.savefig(fig_names[index] + '_avg.png') + else: + fig.savefig(fig_names[index] + '.png') + + +# Annually average the given ROMS data. Note that ROMS takes leap years into +# account for its output, i.e. the 5-day averaging period holds true throughout +# the entire simulation with no days skipped. +# Input: +# time = 1D array of ROMS time values, in years, with year_start added (eg +# 1992.001, 1993.50) assuming year length of 365.25 days +# data = 1D array of ROMS data corresponding to time array +# year_start = integer containing the first year to process +# Output: data_avg = 1D array of annually averaged data +def roms_annual_avg (time, data, year_start): + + data_avg = [] + year = year_start + # Loop over years until we run out of data + while True: + # Find the first index after the beginning of this year + tmp1 = nonzero(time >= year)[0] + if len(tmp1)==0: + # No such index, we have run out of data + break + t_start = tmp1[0] + # Find the first index after the beginning of next year + tmp2 = nonzero(time >= year+1)[0] + if len(tmp2)==0: + # No such index, but we might not have run out of data + # eg simulation that ends on 31 December 2016 + t_end = len(time) + else: + t_end = tmp2[0] + if t_end - t_start < 73: + # This is a partial year, don't count it + break + # Calculate mean between these bounds + data_avg.append(mean(array(data[t_start:t_end]))) + # Increment year + year += 1 + return data_avg + + +# Annually average the given FESOM data. Note that FESOM neglects leap years +# in its output, i.e. the 5-day averaging period will skip 1 day every 4 years, +# so every year has exactly 365/5 = 73 batches of output. +# Input: +# data = 1D array of FESOM data +# days_per_output = averaging period for data +# Output: data_avg = 1D array of annually averaged data +def fesom_annual_avg (data, days_per_output): + + peryear = 365/days_per_output + data_avg = [] + year = 0 + # Loop over years + while True: + if len(data) < peryear*(year+1): + # Run out of data + # FESOM output comes in annual files so we don't have to worry about + # partial years like in ROMS + break + # Calculate mean for this year + data_avg.append(mean(array(data[peryear*year:peryear*(year+1)]))) + # Increment year + year += 1 + return data_avg + + +# Command-line interface +if __name__ == "__main__": + + roms_dir = raw_input("Path to directory containing MetROMS logfiles: ") + fesom_dir = raw_input("Path to directory containing FESOM logfiles: ") + time_flag = int(raw_input("Full timeseries (1) or annual averages (2)? ")) + if time_flag == 1: + annual = False + elif time_flag == 2: + annual = True + mip_timeseries(roms_dir, fesom_dir, annual) + + + + + + + diff --git a/timeseries_massloss.py b/timeseries_massloss.py index 7d86973..c7addc0 100644 --- a/timeseries_massloss.py +++ b/timeseries_massloss.py @@ -61,6 +61,7 @@ def timeseries_massloss (file_path, log_path): # Reached the header for the next ice shelf break index += 1 + f.close() # Calculate dA (masked with ice shelf mask) and lon and lat coordinates print 'Analysing grid'