From 58619399bc1b97cb948bffe86910b6fb93ecc2e4 Mon Sep 17 00:00:00 2001 From: Kaitlin Naughten Date: Mon, 10 Jul 2017 17:42:57 +1000 Subject: [PATCH] Sea ice figures for intercomparison paper --- file_guide.txt | 69 +++-- mip_aice_minmax_nsidc.py | 271 ++++++++++++++---- mip_grid_res.py | 6 +- mip_aice_hi_seasonal.py => mip_hi_seasonal.py | 81 ++---- 4 files changed, 279 insertions(+), 148 deletions(-) rename mip_aice_hi_seasonal.py => mip_hi_seasonal.py (59%) diff --git a/file_guide.txt b/file_guide.txt index d2dac55..90761e9 100644 --- a/file_guide.txt +++ b/file_guide.txt @@ -1234,45 +1234,54 @@ mip_cavity_fields.py: For each major ice shelf, make a 2x1 plot of the given and the FESOM mesh directory. It will output one png for each ice shelf. -mip_aice_hi_seasonal.py: Make two 4x2 plots, showing seasonal averages of - either sea ice concentration or effective thickness - (concentration*height), comparing MetROMS (top) and - FESOM (bottom). - 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 seasonal climatologies of sea ice - concentration and thickness for both MetROMS - and FESOM during the entire simulation, using - seasonal_climatology_cice.py (in roms_tools) - and seasonal_climatology.py (in fesomtools). - Once the files are ready, open python or - ipython and type "run mip_aice_hi_seasonal.py". - The script will prompt you for paths to these - files as well as the FESOM mesh directory. - It will output two png files. - -mip_aice_minmax_nsidc.py: Make a 3x2 plot showing February (top) and August +mip_hi_seasonal.py: Make a 4x2 plot showing seasonal averages of sea ice + effective thickness (concentration*height), comparing + MetROMS (top) and FESOM (bottom). + 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 seasonal + climatologies of sea ice concentration and + thickness for both MetROMS and FESOM during the + entire simulation, using + seasonal_climatology_cice.py (in roms_tools) and + seasonal_climatology.py (in fesomtools). Once the + files are ready, open python or ipython and type + "run mip_aice_hi_seasonal.py". The script will + prompt you for paths to these files as well as the + FESOM mesh directory. It will output a png file. + +mip_aice_minmax_nsidc.py: Make a 4x2 plot showing February (top) and September (bottom) sea ice concentration (1992-2015 average) for - NSIDC, MetROMS, and FESOM. + NSIDC, MetROMS, and FESOM. The fourth column shows + timeseries of February and September sea ice extent + for NSIDC, MetROMS, and FESOM. 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, + the cloned repository on your system. Make + sure timeseries_seaice_exent.py, in both + roms_tools and fesomtools, has been run for + the MetROMS and FESOM simulations. Next, download the NSIDC Climate Data Record monthly averages of Southern Hemisphere sea ice - concentration, and make sure the variables - nsidc_* near the top of the file will - reconstruct the right filenames on your - system. Then open python or ipython and type + concentration, and the NSIDC Sea Ice Index + csv files with total Antarctic sea ice extent + for each year, in February and September. + Make sure the variables nsidc_* near the top + of the file will reconstruct the right + filenames on your system. Then open python or + ipython and type "run mip_aice_minmax_nsidc.py". The script will prompt you for paths to the CICE file containing 5-day averages for the entire - simulation, the FESOM mesh directory, and the - FESOM output directory containing one - ice.mean.nc file (5-day averages) for each - year. It will output a png file. + simulation, the CICE logfile for + timeseries_seaice_extent.py, the FESOM mesh + directory, the FESOM output directory + containing one ice.mean.nc file (5-day + averages) for each year, and the FESOM + logfile for timeseries_seaice_extent.py in + fesomtools. It will output a png file. mip_zonal_cavity_ts.py: Create one 2x2 plot for each major ice shelf: zonal slices of temperature (top) and salinity (bottom), diff --git a/mip_aice_minmax_nsidc.py b/mip_aice_minmax_nsidc.py index cfe4a66..c6f895e 100644 --- a/mip_aice_minmax_nsidc.py +++ b/mip_aice_minmax_nsidc.py @@ -1,4 +1,4 @@ -from netCDF4 import Dataset +from netCDF4 import Dataset, num2date from numpy import * from matplotlib.collections import PatchCollection from matplotlib.pyplot import * @@ -9,19 +9,28 @@ from patches import * from monthly_avg import * -# Make a 3x2 plot showing February (top) and August (bottom) sea ice -# concentration (1992-2015 average) for NSIDC, MetROMS, and FESOM. +# Make a 4x2 plot showing February (top) and September (bottom) sea ice +# concentration (1992-2015 average) for NSIDC, MetROMS, and FESOM. The fourth +# column shows timeseries of February and September sea ice extent for NSIDC, +# MetROMS, and FESOM. # Input: # cice_file = path to CICE output file containing 5-day averages for the entire # simulation +# cice_log = path to CICE logfile from timeseries_seaice_extent.py # fesom_mesh_path = path to FESOM mesh directory # fesom_output_dir = path to FESOM output directory containing one ice.mean.nc # file for each year (5-day averages) -def mip_aice_minmax_nsidc (cice_file, fesom_mesh_path, fesom_output_dir): +# fesom_log = path to FESOM logfile from timeseries_seaice_extent.py in +# fesomtools +def mip_aice_minmax_nsidc (cice_file, cice_log, fesom_mesh_path, fesom_output_dir, fesom_log): # Range of years to process (exclude 2016 because no NSIDC data) start_year = 1992 end_year = 2015 + # Starting and ending days for each month + # Ignore leap years, they will be dealt with later + start_day = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] + end_day = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] # Beginning of FESOM output filenames fesom_file_head = 'MK44005' # Paths to reconstruct NSIDC files for each month @@ -29,6 +38,9 @@ def mip_aice_minmax_nsidc (cice_file, fesom_mesh_path, fesom_output_dir): nsidc_head2 = '/short/m68/kaa561/nsidc_aice/seaice_conc_monthly_sh_f13_' nsidc_head3 = '/short/m68/kaa561/nsidc_aice/seaice_conc_monthly_sh_f17_' nsidc_tail = '_v02r00.nc' + # Path to NSIDC monthly extent timeseries for February and September + nsidc_feb_ts_file = '/short/m68/kaa561/nsidc_aice/S_02_extent_v2.1.csv' + nsidc_sep_ts_file = '/short/m68/kaa561/nsidc_aice/S_09_extent_v2.1.csv' # FESOM plotting parameters circumpolar = True mask_cavities = True @@ -37,26 +49,27 @@ def mip_aice_minmax_nsidc (cice_file, fesom_mesh_path, fesom_output_dir): num_years = end_year - start_year + 1 print 'Processing NSIDC' + # First read the grid id = Dataset(nsidc_head1 + '199201' + nsidc_tail, 'r') nsidc_lon = id.variables['longitude'][:,:] nsidc_lat = id.variables['latitude'][:,:] id.close() - # Read February and August concentration for each year + # Read February and September concentration for each year nsidc_feb = ma.empty([num_years, size(nsidc_lon,0), size(nsidc_lat,1)]) - nsidc_aug = ma.empty([num_years, size(nsidc_lon,0), size(nsidc_lat,1)]) + nsidc_sep = ma.empty([num_years, size(nsidc_lon,0), size(nsidc_lat,1)]) # Loop over years for year in range(start_year, end_year+1): # Reconstruct file paths if year < 1996: feb_file = nsidc_head1 + str(year) + '02' + nsidc_tail - aug_file = nsidc_head1 + str(year) + '08' + nsidc_tail + sep_file = nsidc_head1 + str(year) + '09' + nsidc_tail elif year < 2008: feb_file = nsidc_head2 + str(year) + '02' + nsidc_tail - aug_file = nsidc_head2 + str(year) + '08' + nsidc_tail + sep_file = nsidc_head2 + str(year) + '09' + nsidc_tail else: feb_file = nsidc_head3 + str(year) + '02' + nsidc_tail - aug_file = nsidc_head3 + str(year) + '08' + nsidc_tail + sep_file = nsidc_head3 + str(year) + '09' + nsidc_tail # Read February data and mask id = Dataset(feb_file, 'r') feb_data_tmp = id.variables['seaice_conc_monthly_cdr'][0,:,:] @@ -70,22 +83,22 @@ def mip_aice_minmax_nsidc (cice_file, fesom_mesh_path, fesom_output_dir): feb_data[nsidc_mask.mask] = ma.masked # Save result nsidc_feb[year-start_year,:,:] = feb_data[:,:] - # Repeat for August - id = Dataset(aug_file, 'r') - aug_data_tmp = id.variables['seaice_conc_monthly_cdr'][0,:,:] + # Repeat for September + id = Dataset(sep_file, 'r') + sep_data_tmp = id.variables['seaice_conc_monthly_cdr'][0,:,:] nsidc_mask = id.variables['stdev_of_seaice_conc_monthly_cdr'][0,:,:] id.close() - aug_data = ma.empty(shape(aug_data_tmp)) - aug_data[:,:] = 0.0 - aug_data[~nsidc_mask.mask] = aug_data_tmp[~nsidc_mask.mask] - aug_data[nsidc_mask.mask] = ma.masked - nsidc_aug[year-start_year,:,:] = aug_data[:,:] + sep_data = ma.empty(shape(sep_data_tmp)) + sep_data[:,:] = 0.0 + sep_data[~nsidc_mask.mask] = sep_data_tmp[~nsidc_mask.mask] + sep_data[nsidc_mask.mask] = ma.masked + nsidc_sep[year-start_year,:,:] = sep_data[:,:] # Average over years nsidc_feb = mean(nsidc_feb, axis=0) - nsidc_aug = mean(nsidc_aug, axis=0) + nsidc_sep = mean(nsidc_sep, axis=0) # Make sure mask is still there nsidc_feb[nsidc_mask.mask] = ma.masked - nsidc_aug[nsidc_mask.mask] = ma.masked + nsidc_sep[nsidc_mask.mask] = ma.masked # Polar coordinates for plotting nsidc_x = -(nsidc_lat+90)*cos(nsidc_lon*deg2rad+pi/2) nsidc_y = (nsidc_lat+90)*sin(nsidc_lon*deg2rad+pi/2) @@ -95,7 +108,32 @@ def mip_aice_minmax_nsidc (cice_file, fesom_mesh_path, fesom_output_dir): bdry3 = amin(nsidc_y[:,0]) bdry4 = amax(nsidc_y[:,-1]) + # Now read extent timeseries + nsidc_feb_extent = [] + f = open(nsidc_feb_ts_file, 'r') + f.readline() + # Skip the years we don't care about + for year in range(1979, start_year): + f.readline() + # Read the years we care about + for year in range(start_year, end_year+1): + tmp = f.readline().split(',') + # Extract the extent (second last column) + nsidc_feb_extent.append(float(tmp[-2])) + f.close() + # Repeat for September + nsidc_sep_extent = [] + f = open(nsidc_sep_ts_file, 'r') + f.readline() + for year in range(1979, start_year): + f.readline() + for year in range(start_year, end_year+1): + tmp = f.readline().split(',') + nsidc_sep_extent.append(float(tmp[-2])) + f.close() + print 'Processing MetROMS' + # First read the grid id = Dataset(cice_file, 'r') cice_lon_tmp = id.variables['TLON'][:,:] @@ -108,61 +146,163 @@ def mip_aice_minmax_nsidc (cice_file, fesom_mesh_path, fesom_output_dir): cice_lon[:,-1] = cice_lon_tmp[:,0] cice_lat[:,:-1] = cice_lat_tmp cice_lat[:,-1] = cice_lat_tmp[:,0] - # Get averages for February and August + # Get averages for February and September # Start with first year just to initialise the arrays with the right size print '...monthly average for ' + str(start_year) cice_feb_tmp = monthly_avg_cice(cice_file, 'aice', shape(cice_lon_tmp), 1, instance=1) - cice_aug_tmp = monthly_avg_cice(cice_file, 'aice', shape(cice_lon_tmp), 7, instance=1) + cice_sep_tmp = monthly_avg_cice(cice_file, 'aice', shape(cice_lon_tmp), 8, instance=1) # Loop over the rest of the years for year in range(start_year+1, end_year+1): print '...monthly average for ' + str(year) cice_feb_tmp = cice_feb_tmp + monthly_avg_cice(cice_file, 'aice', shape(cice_lon_tmp), 1, instance=year-start_year+1) - cice_aug_tmp = cice_aug_tmp + monthly_avg_cice(cice_file, 'aice', shape(cice_lon_tmp), 7, instance=year-start_year+1) + cice_sep_tmp = cice_sep_tmp + monthly_avg_cice(cice_file, 'aice', shape(cice_lon_tmp), 8, instance=year-start_year+1) # Convert from integrals to averages cice_feb_tmp = cice_feb_tmp/num_years - cice_aug_tmp = cice_aug_tmp/num_years + cice_sep_tmp = cice_sep_tmp/num_years # Wrap the periodic boundary cice_feb = ma.empty(shape(cice_lon)) - cice_aug = ma.empty(shape(cice_lon)) + cice_sep = ma.empty(shape(cice_lon)) cice_feb[:,:-1] = cice_feb_tmp - cice_aug[:,:-1] = cice_aug_tmp + cice_sep[:,:-1] = cice_sep_tmp cice_feb[:,-1] = cice_feb_tmp[:,0] - cice_aug[:,-1] = cice_aug_tmp[:,0] + cice_sep[:,-1] = cice_sep_tmp[:,0] # Polar coordinates for plotting cice_x = -(cice_lat+90)*cos(cice_lon*deg2rad+pi/2) cice_y = (cice_lat+90)*sin(cice_lon*deg2rad+pi/2) + # Now get extent timeseries + # Read 5-day logfile + cice_time_vals = [] + cice_extent_5day = [] + f = open(cice_log, 'r') + f.readline() + for line in f: + try: + cice_time_vals.append(float(line)) + except(ValueError): + break + for line in f: + cice_extent_5day.append(float(line)) + f.close() + # Convert time to Date objects + cice_time = num2date(array(cice_time_vals)*365.25, units='days since 1992-01-01 00:00:00', calendar='gregorian') + # Initialise integral arrays for monthly averages + # Add an extra year because the simulation goes to the end of 2016 + cice_extent = zeros((num_years+1)*12) + cice_ndays = zeros((num_years+1)*12) + for t in range(size(cice_time)): + # 5-day averages marked with the next day's date + year = cice_time[t].year + month = cice_time[t].month-1 # Convert to 0-indexed + day = cice_time[t].day + # Check for leap years + leap_year = False + if mod(year, 4) == 0: + leap_year = True + if mod(year, 100) == 0: + leap_year = False + if mod(year, 400) == 0: + leap_year = True + if leap_year: + end_day[1] = 29 + else: + end_day[1] = 28 + if day-5 < start_day[month]: + # Spills over into the previous month + prev_month = mod(month-1, 12) + # How many days does it spill over by? + spill_days = start_day[month]-day+5 + # Should be between 1 and 5 + if spill_days < 1 or spill_days > 5: + print 'Problem: spill_days is ' + str(spill_days) + print 'Timestep ' + str(t+1) + print 'Year ' + str(year) + print 'Month ' + str(month+1) + print 'Day ' + str(day) + #return + # Split between previous month and this month + # First find indices to update + if prev_month+1 == 12: + # Spilled into previous year + index_prev = (year-1-start_year)*12 + prev_month + else: + index_prev = (year-start_year)*12 + prev_month + index = (year-start_year)*12 + month + # Integrate + cice_extent[index_prev] += cice_extent_5day[t]*spill_days + cice_ndays[index_prev] += spill_days + cice_extent[index] += cice_extent_5day[t]*(5-spill_days) + cice_ndays[index] += 5-spill_days + else: + # Entirely within the month + index = (year-start_year)*12 + month + cice_extent[index] += cice_extent_5day[t]*5 + cice_ndays[index] += 5 + # Convert from integrals to averages + cice_extent /= cice_ndays + # Extract February and September + cice_feb_extent = [] + cice_sep_extent = [] + for year in range(start_year, end_year+1): + cice_feb_extent.append(cice_extent[(year-start_year)*12+1]) + cice_sep_extent.append(cice_extent[(year-start_year)*12+8]) + print 'Processing FESOM' + # First build the grid elements, patches = make_patches(fesom_mesh_path, circumpolar, mask_cavities) - # Get averages for February and August + # Get averages for February and September # Start with first year just to initialise the arrays with the right size print '...monthly average for ' + str(start_year) fesom_feb_nodes = monthly_avg(fesom_output_dir + fesom_file_head + '.' + str(start_year) + '.ice.mean.nc', 'area', 1) - fesom_aug_nodes = monthly_avg(fesom_output_dir + fesom_file_head + '.' + str(start_year) + '.ice.mean.nc', 'area', 7) + fesom_sep_nodes = monthly_avg(fesom_output_dir + fesom_file_head + '.' + str(start_year) + '.ice.mean.nc', 'area', 8) # Loop over the rest of the years for year in range(start_year+1, end_year+1): print '...monthly average for ' + str(year) fesom_feb_nodes = fesom_feb_nodes + monthly_avg(fesom_output_dir + fesom_file_head + '.' + str(year) + '.ice.mean.nc', 'area', 1) - fesom_aug_nodes = fesom_aug_nodes + monthly_avg(fesom_output_dir + fesom_file_head + '.' + str(year) + '.ice.mean.nc', 'area', 7) + fesom_sep_nodes = fesom_sep_nodes + monthly_avg(fesom_output_dir + fesom_file_head + '.' + str(year) + '.ice.mean.nc', 'area', 8) # Convert from integrals to averages fesom_feb_nodes = fesom_feb_nodes/num_years - fesom_aug_nodes = fesom_aug_nodes/num_years + fesom_sep_nodes = fesom_sep_nodes/num_years # Find element-averages fesom_feb = [] - fesom_aug = [] + fesom_sep = [] for elm in elements: if not elm.cavity: # Average over 3 component nodes fesom_feb.append(mean(array([fesom_feb_nodes[elm.nodes[0].id], fesom_feb_nodes[elm.nodes[1].id], fesom_feb_nodes[elm.nodes[2].id]]))) - fesom_aug.append(mean(array([fesom_aug_nodes[elm.nodes[0].id], fesom_aug_nodes[elm.nodes[1].id], fesom_aug_nodes[elm.nodes[2].id]]))) + fesom_sep.append(mean(array([fesom_sep_nodes[elm.nodes[0].id], fesom_sep_nodes[elm.nodes[1].id], fesom_sep_nodes[elm.nodes[2].id]]))) fesom_feb = array(fesom_feb) - fesom_aug = array(fesom_aug) + fesom_sep = array(fesom_sep) + + # Get extent timeseries + # Read 5-day logfile + fesom_extent_5day = [] + f = open(fesom_log, 'r') + f.readline() + for line in f: + fesom_extent_5day.append(float(line)) + f.close() + # Initialise monthly arrays + fesom_feb_extent = zeros(num_years) + fesom_sep_extent = zeros(num_years) + for year in range(start_year, end_year+1): + # First timestep of year in 5-day logfile + t0 = (year-start_year)*73 + # Monthly averages are hard-coded and ugly + # Feburary: 4/5 of index 7, indices 8-11, and 4/5 of index 12 + fesom_feb_extent[year-start_year] = (fesom_extent_5day[t0+6]*4 + sum(fesom_extent_5day[t0+7:t0+11]*5) + fesom_extent_5day[t0+11]*4)/28.0 + # September: 2/5 of index 49, indices 50-54, 3/5 of index 55 + fesom_sep_extent[year-start_year] = (fesom_extent_5day[t0+48]*2 + sum(fesom_extent_5day[t0+49:t0+54]*5) + fesom_extent_5day[t0+54]*3)/30.0 + + time_axis = arange(start_year, end_year+1) print 'Plotting' - fig = figure(figsize=(16,10)) + fig = figure(figsize=(20,10)) + gs1 = GridSpec(2, 3) + gs1.update(left=0.12, right=0.7, wspace=0.05, hspace=0.05) # NSIDC, February - ax = fig.add_subplot(2, 3, 1, aspect='equal') + ax = subplot(gs1[0, 0], aspect='equal') img = pcolor(nsidc_x, nsidc_y, nsidc_feb, vmin=0, vmax=1, cmap='jet') xlim([bdry1, bdry2]) ylim([bdry3, bdry4]) @@ -170,14 +310,14 @@ def mip_aice_minmax_nsidc (cice_file, fesom_mesh_path, fesom_output_dir): title('NSIDC', fontsize=24) text(-39, 0, 'February', fontsize=24, ha='right') # MetROMS, February - ax = fig.add_subplot(2, 3, 2, aspect='equal') + ax = subplot(gs1[0, 1], aspect='equal') img = pcolor(cice_x, cice_y, cice_feb, vmin=0, vmax=1, cmap='jet') xlim([bdry1, bdry2]) ylim([bdry3, bdry4]) axis('off') title('MetROMS', fontsize=24) # FESOM, February - ax = fig.add_subplot(2, 3, 3, aspect='equal') + ax = subplot(gs1[0, 2], aspect='equal') img = PatchCollection(patches, cmap='jet') img.set_array(fesom_feb) img.set_clim(vmin=0, vmax=1) @@ -186,24 +326,26 @@ def mip_aice_minmax_nsidc (cice_file, fesom_mesh_path, fesom_output_dir): xlim([bdry1, bdry2]) ylim([bdry3, bdry4]) axis('off') - title('FESOM', fontsize=24) - # NSIDC, August - ax = fig.add_subplot(2, 3, 4, aspect='equal') - img = pcolor(nsidc_x, nsidc_y, nsidc_aug, vmin=0, vmax=1, cmap='jet') + title('FESOM (high-res)', fontsize=24) + # Main title + text(-170, 47, 'a) Sea ice concentration ('+str(start_year)+'-'+str(end_year)+' average)', fontsize=30) + # NSIDC, September + ax = subplot(gs1[1, 0], aspect='equal') + img = pcolor(nsidc_x, nsidc_y, nsidc_sep, vmin=0, vmax=1, cmap='jet') xlim([bdry1, bdry2]) ylim([bdry3, bdry4]) axis('off') - text(-39, 0, 'August', fontsize=24, ha='right') - # MetROMS, August - ax = fig.add_subplot(2, 3, 5, aspect='equal') - img = pcolor(cice_x, cice_y, cice_aug, vmin=0, vmax=1, cmap='jet') + text(-39, 0, 'September', fontsize=24, ha='right') + # MetROMS, September + ax = subplot(gs1[1, 1], aspect='equal') + img = pcolor(cice_x, cice_y, cice_sep, vmin=0, vmax=1, cmap='jet') xlim([bdry1, bdry2]) ylim([bdry3, bdry4]) axis('off') - # FESOM, August - ax = fig.add_subplot(2, 3, 6, aspect='equal') + # FESOM, September + ax = subplot(gs1[1, 2], aspect='equal') img = PatchCollection(patches, cmap='jet') - img.set_array(fesom_aug) + img.set_array(fesom_sep) img.set_clim(vmin=0, vmax=1) img.set_edgecolor('face') ax.add_collection(img) @@ -211,13 +353,32 @@ def mip_aice_minmax_nsidc (cice_file, fesom_mesh_path, fesom_output_dir): ylim([bdry3, bdry4]) axis('off') # Add a colourbar at the bottom - cbaxes = fig.add_axes([0.35, 0.04, 0.3, 0.04]) + cbaxes = fig.add_axes([0.26, 0.04, 0.3, 0.04]) cbar = colorbar(img, orientation='horizontal', ticks=arange(0,1+0.25,0.25), cax=cbaxes) cbar.ax.tick_params(labelsize=16) - # Main title - suptitle('Sea ice concentration ('+str(start_year)+'-'+str(end_year)+' average)', fontsize=30) - # Make panels closer together - subplots_adjust(wspace=0.05, hspace=0.05) + # Add extent timeseries on rightmost column, with more space for labels + gs2 = GridSpec(2, 1) + gs2.update(left=0.73, right=0.95, wspace=0.1, hspace=0.15) + # February + ax = subplot(gs2[0, 0]) + ax.plot(time_axis, nsidc_feb_extent, color='red', linewidth=1.5) + ax.plot(time_axis, cice_feb_extent, color='blue', linewidth=1.5) + ax.plot(time_axis, fesom_feb_extent, color='green', linewidth=1.5) + xlim([start_year, end_year]) + ax.tick_params(axis='x', labelsize=16) + ax.tick_params(axis='y', labelsize=16) + grid(True) + title('b) Sea ice extent\n'+r'(million km$^2$)', fontsize=26) + # Extent timeseries, September + ax = subplot(gs2[1, 0]) + ax.plot(time_axis, nsidc_sep_extent, color='red', label='NSIDC', linewidth=1.5) + ax.plot(time_axis, cice_sep_extent, color='blue', label='MetROMS', linewidth=1.5) + ax.plot(time_axis, fesom_sep_extent, color='green', label='FESOM (high-res)', linewidth=1.5) + xlim([start_year, end_year]) + ax.tick_params(axis='x', labelsize=16) + ax.tick_params(axis='y', labelsize=16) + grid(True) + ax.legend(bbox_to_anchor=(1.15,-0.1), ncol=3) #fig.show() fig.savefig('aice_minmax_nsidc.png') @@ -226,9 +387,11 @@ def mip_aice_minmax_nsidc (cice_file, fesom_mesh_path, fesom_output_dir): if __name__ == "__main__": cice_file = raw_input("Path to CICE output file containing data for the entire simulation: ") + cice_log = raw_input("Path to CICE sea ice extent logfile: ") fesom_mesh_path = raw_input("Path to FESOM mesh directory: ") fesom_output_dir = raw_input("Path to FESOM output directory containing one ice.mean.nc file for each year: ") - mip_aice_minmax_nsidc (cice_file, fesom_mesh_path, fesom_output_dir) + fesom_log = raw_input("Path to FESOM sea ice extent logfile: ") + mip_aice_minmax_nsidc(cice_file, cice_log, fesom_mesh_path, fesom_output_dir, fesom_log) diff --git a/mip_grid_res.py b/mip_grid_res.py index 7adccf5..5dfc7e2 100644 --- a/mip_grid_res.py +++ b/mip_grid_res.py @@ -70,7 +70,7 @@ def mip_grid_res (roms_grid_file, fesom_mesh_low, fesom_mesh_high, save=False, f xlim([-lat_max, lat_max]) ylim([-lat_max, lat_max]) axis('off') - title('MetROMS', fontsize=28) + title('a) MetROMS', fontsize=28) # FESOM low-res ax2 = fig.add_subplot(1,3,2, aspect='equal') img_low = PatchCollection(patches_low, cmap='jet') @@ -81,7 +81,7 @@ def mip_grid_res (roms_grid_file, fesom_mesh_low, fesom_mesh_high, save=False, f xlim([-lat_max, lat_max]) ylim([-lat_max, lat_max]) axis('off') - title('FESOM low-res', fontsize=28) + title('b) FESOM low-res', fontsize=28) # FESOM high-res ax3 = fig.add_subplot(1,3,3, aspect='equal') img_high = PatchCollection(patches_high, cmap='jet') @@ -92,7 +92,7 @@ def mip_grid_res (roms_grid_file, fesom_mesh_low, fesom_mesh_high, save=False, f xlim([-lat_max, lat_max]) ylim([-lat_max, lat_max]) axis('off') - title('FESOM high-res', fontsize=28) + title('c) FESOM high-res', fontsize=28) cbaxes = fig.add_axes([0.92, 0.2, 0.01, 0.6]) cbar = colorbar(img_high, cax=cbaxes, extend='max', ticks=arange(limits[0], limits[1]+5, 5)) cbar.ax.tick_params(labelsize=24) diff --git a/mip_aice_hi_seasonal.py b/mip_hi_seasonal.py similarity index 59% rename from mip_aice_hi_seasonal.py rename to mip_hi_seasonal.py index 4228023..6d466a9 100644 --- a/mip_aice_hi_seasonal.py +++ b/mip_hi_seasonal.py @@ -7,17 +7,16 @@ sys.path.insert(0, '/short/y99/kaa561/fesomtools') from patches import * -# Make two 4x2 plots, showing seasonal averages of either sea ice concentration -# or effective thickness (concentration*height), comparing MetROMS (top) and -# FESOM (bottom). +# Make a 4x2 plot showing seasonal averages of sea ice effective thickness +# (concentration*height), comparing MetROMS (top) and FESOM (bottom). # Input: # cice_seasonal_file = path to seasonal climatology of CICE variables 'aice' and # 'hi', pre-computed using seasonal_climatology_cice.py -# fesom_mesh_path = path to FESOM mesh directory -# fesom_seasonal_file = path to seasonal climatology of FESOM variables 'area' -# and 'hice', pre-computed using seasonal_climatology.py -# in the "fesomtools" repository -def mip_aice_hi_seasonal (cice_seasonal_file, fesom_mesh_path, fesom_seasonal_file): +# fesom_mesh_path = path to FESOM mesh directories +# fesom_seasonal_file = path to seasonal climatology of FESOM variable 'hice', +# pre-computed using seasonal_climatology.py in the +# "fesomtools" repository +def mip_hi_seasonal (cice_seasonal_file, fesom_mesh_path, fesom_seasonal_file): # Northern boundary of plot 50S nbdry = -50 + 90 @@ -29,8 +28,7 @@ def mip_aice_hi_seasonal (cice_seasonal_file, fesom_mesh_path, fesom_seasonal_fi # Season names for plot titles season_names = ['DJF', 'MAM', 'JJA', 'SON'] # Colour bounds - aice_bounds = [0, 1] - hi_bounds = [0, 1.5] + bounds = [0, 1.5] print 'Processing MetROMS' # Read CICE grid data @@ -65,7 +63,6 @@ def mip_aice_hi_seasonal (cice_seasonal_file, fesom_mesh_path, fesom_seasonal_fi elements, patches = make_patches(fesom_mesh_path, circumpolar, mask_cavities) # Read data id = Dataset(fesom_seasonal_file, 'r') - fesom_aice_nodes = id.variables['area'][:,:] fesom_hi_nodes = id.variables['hice'][:,:] id.close() # Count the number of elements not in ice shelf cavities @@ -73,59 +70,22 @@ def mip_aice_hi_seasonal (cice_seasonal_file, fesom_mesh_path, fesom_seasonal_fi for elm in elements: if not elm.cavity: num_elm += 1 - # Set up arrays for element-averages for each season - fesom_aice = zeros([4, num_elm]) + # Set up array for element-averages for each season fesom_hi = zeros([4, num_elm]) - # Loop over elements to fill these in + # Loop over elements to fill this in i = 0 for elm in elements: if not elm.cavity: # Average over 3 component nodes - fesom_aice[:,i] = (fesom_aice_nodes[:,elm.nodes[0].id] + fesom_aice_nodes[:,elm.nodes[1].id] + fesom_aice_nodes[:,elm.nodes[2].id])/3 fesom_hi[:,i] = (fesom_hi_nodes[:,elm.nodes[0].id] + fesom_hi_nodes[:,elm.nodes[1].id] + fesom_hi_nodes[:,elm.nodes[2].id])/3 i += 1 - print 'Plotting sea ice concentration' - fig = figure(figsize=(20,9)) - # Loop over seasons + print 'Plotting' + fig = figure(figsize=(19,9)) for season in range(4): # MetROMS ax = fig.add_subplot(2, 4, season+1, aspect='equal') - pcolor(cice_x, cice_y, cice_aice[season,:,:], vmin=aice_bounds[0], vmax=aice_bounds[1], cmap='jet') - if season == 0: - text(-43, 0, 'MetROMS', fontsize=24, ha='right') - title(season_names[season], fontsize=24) - xlim([-nbdry, nbdry]) - ylim([-nbdry, nbdry]) - axis('off') - # FESOM - ax = fig.add_subplot(2, 4, season+5, aspect='equal') - img = PatchCollection(patches, cmap='jet') - img.set_array(fesom_aice[season,:]) - img.set_clim(vmin=aice_bounds[0], vmax=aice_bounds[1]) - img.set_edgecolor('face') - ax.add_collection(img) - xlim([-nbdry, nbdry]) - ylim([-nbdry, nbdry]) - axis('off') - if season == 0: - text(-43, 0, 'FESOM', fontsize=24, ha='right') - # Add a horizontal colorbar at the bottom - cbaxes = fig.add_axes([0.25, 0.04, 0.5, 0.02]) - cbar = colorbar(img, orientation='horizontal', ticks=arange(aice_bounds[0],aice_bounds[1]+0.25,0.25), cax=cbaxes) - cbar.ax.tick_params(labelsize=16) - # Add the main title - suptitle('Sea ice concentration', fontsize=30) - # Decrease space between plots - subplots_adjust(wspace=0.025,hspace=0.025) - fig.savefig('aice_seasonal.png') - - print 'Plotting sea ice thickness' - fig = figure(figsize=(20,9)) - for season in range(4): - # MetROMS - ax = fig.add_subplot(2, 4, season+1, aspect='equal') - pcolor(cice_x, cice_y, cice_hi[season,:,:], vmin=hi_bounds[0], vmax=hi_bounds[1], cmap='jet') + pcolor(cice_x, cice_y, cice_hi[season,:,:], vmin=bounds[0], vmax=bounds[1], cmap='jet') if season == 0: text(-43, 0, 'MetROMS', fontsize=24, ha='right') title(season_names[season], fontsize=24) @@ -136,7 +96,7 @@ def mip_aice_hi_seasonal (cice_seasonal_file, fesom_mesh_path, fesom_seasonal_fi ax = fig.add_subplot(2, 4, season+5, aspect='equal') img = PatchCollection(patches, cmap='jet') img.set_array(fesom_hi[season,:]) - img.set_clim(vmin=hi_bounds[0], vmax=hi_bounds[1]) + img.set_clim(vmin=bounds[0], vmax=bounds[1]) img.set_edgecolor('face') ax.add_collection(img) xlim([-nbdry, nbdry]) @@ -144,9 +104,10 @@ def mip_aice_hi_seasonal (cice_seasonal_file, fesom_mesh_path, fesom_seasonal_fi axis('off') if season == 0: text(-43, 0, 'FESOM', fontsize=24, ha='right') - cbaxes = fig.add_axes([0.25, 0.04, 0.5, 0.02]) - cbar = colorbar(img, orientation='horizontal', ticks=arange(hi_bounds[0],hi_bounds[1]+0.5,0.5), cax=cbaxes, extend='max') - cbar.ax.tick_params(labelsize=16) + text(-43, -10, '(high-res)', fontsize=24,ha='right') + cbaxes = fig.add_axes([0.35, 0.04, 0.3, 0.02]) + cbar = colorbar(img, orientation='horizontal', ticks=arange(bounds[0],bounds[1]+0.5,0.5), cax=cbaxes, extend='max') + cbar.ax.tick_params(labelsize=20) suptitle('Sea ice effective thickness (m)', fontsize=30) subplots_adjust(wspace=0.025,hspace=0.025) fig.savefig('hi_seasonal.png') @@ -157,8 +118,6 @@ def mip_aice_hi_seasonal (cice_seasonal_file, fesom_mesh_path, fesom_seasonal_fi cice_seasonal_file = raw_input("Path to CICE seasonal climatology file containing aice and hi: ") fesom_mesh_path = raw_input("Path to FESOM mesh directory: ") - fesom_seasonal_file = raw_input("Path to FESOM seasonal climatology file containing area and hice: ") - mip_aice_hi_seasonal(cice_seasonal_file, fesom_mesh_path, fesom_seasonal_file) - - + fesom_seasonal_file = raw_input("Path to FESOM seasonal climatology file containing hice: ") + mip_hi_seasonal(cice_seasonal_file, fesom_mesh_path, fesom_seasonal_file)