Added plotting scripts for mass loss error map and NSIDC sea ice comp…

…arison

README.md

@@ -1,4 +1,4 @@
-Python scripts for pre- and post-processing of ROMS simulations. See file_guide.txt for descriptions of what they do and how to use them.
+Python scripts for pre- and post-processing of ROMS simulations, and some interactive plotting routines. See file_guide.txt for descriptions of what they do and how to use them.
 
 Created for the circumpolar Antarctic configuration of ROMS-CICE-MCT with a quarter-degree grid; should be easy to adapt to other domains.
 

aice_animation.py

@@ -4,8 +4,7 @@
 from matplotlib.animation import *
 
 # Create an animation of sea ice concentration for the given simulation.
-# Save as an mp4 file. (This doesn't run well on the CCRC desktops so use
-# something like a MacBook instead.)
+# Save as an mp4 file.
 # In order for the mp4 saving to work, you must first type "module load ffmpeg"
 # on raijin before opening ipython.
 # Couldn't work out how to make this an encapsulated function so it is just
@@ -14,9 +13,9 @@
 # Directory containing CICE output files
 directory = '/short/y99/kaa561/roms_spinup_newest/cice/'
 # Number of time indices in each file
-num_ts = [288]
+num_ts = [288, 288, 270, 288, 270, 90]
 # File number to start with for the animation (1-based)
-start_file = 1
+start_file = 6
 # Degrees to radians conversion factor
 deg2rad = pi/180
 # Names of each month for making titles
@@ -66,6 +65,6 @@ def animate(i):
     return img
 
 # Animate once every time index from start_file to the last file
-anim = FuncAnimation(fig, func=animate, frames=range(215,288)) #sum(array(num_ts[start_file-1:])))
+anim = FuncAnimation(fig, func=animate, frames=range(17,90)) #sum(array(num_ts[start_file-1:])))
 # Save as an mp4 with one frame per second
 anim.save('aice.mp4', fps=1)

cmip5_eraint_rms_errors.py

@@ -22,7 +22,7 @@ def cmip5_eraint_rms_errors():
     # Variable names in NetCDF files
     var_names = ['Pair', 'Tair', 'Hair', 'cloud', 'Uwind', 'Vwind', 'precip', 'snow', 'evap', 'swrad', 'lwrad']
     # Path to ROMS grid
-    roms_grid = '/short/m68/kaa561/ROMS-CICE-MCT/apps/common/grid/circ30S_quarterdegree_rp5.nc'
+    roms_grid = '/short/m68/kaa561/ROMS-CICE-MCT/apps/common/grid/circ30S_quarterdegree_10m.nc'
     # Radius of the Earth in metres
     r = 6.371e6
     # Degrees to radians conversion factor

cmip5_plot.py

@@ -31,7 +31,7 @@ def cmip5_plot (var, season, model_names, save=False, fig_name=None):
     # and "ocean/"
     directory = '/short/y99/kaa561/CMIP5_forcing/'
     # Path to ROMS grid file
-    roms_grid = '/short/m68/kaa561/ROMS-CICE-MCT/apps/common/grid/circ30S_quarterdegree_rp5.nc'
+    roms_grid = '/short/m68/kaa561/ROMS-CICE-MCT/apps/common/grid/circ30S_quarterdegree_10m.nc'
 
     # Figure out whether this is an atmosphere or ocean variable
     if var in ['Pair', 'Tair', 'Hair', 'cloud', 'Uwind', 'Vwind', 'precip', 'snow', 'evap', 'swrad', 'lwrad']:

convert_era.job

@@ -1,8 +1,8 @@
 #!/bin/bash
 #PBS -N convert_era
-#PBS -P gh8
-#PBS -q normal
-#PBS -l walltime=3:00:00,ncpus=1,mem=31gb
+#PBS -P y99
+#PBS -q express
+#PBS -l walltime=3:00:00,ncpus=1,mem=30gb
 #PBS -j oe
 #PBS -v YEAR,COUNT
 
@@ -13,6 +13,11 @@
 # To get it started for eg 1992, type
 # qsub -v YEAR=1992,COUNT=0 convert_era.job
 
+module unload python/2.7.3
+module unload python/2.7.3-matplotlib
+module load python/2.7.6
+module load python/2.7.6-matplotlib
+
 echo "YEAR = $YEAR"
 echo "COUNT = $COUNT"
 cd $PBS_O_WORKDIR

file_guide.txt

@@ -162,11 +162,12 @@ massloss.py: Calculates and plots timeseries of basal mass loss and
 	     have to be recomputed if the run is extended; at the beginning of
 	     this script, previous values will be read from the same log file
 	     if it exists.
-	     To run: First edit the limits on i- and j- coordinates for each
-	             ice shelf to suit your grid. Then open python or ipython,
-		     and type "run massloss.py". The script will prompt you for
-		     the paths to the ROMS grid file, the ocean history or
-		     averages file, and the log file.
+	     To run: Open python or ipython, and type "run massloss.py". The
+	             script will prompt you for the paths to the ROMS grid file,
+		     the ocean history or averages file, and the log file. If
+		     you are using ice shelf draft data from something other
+		     than RTopo 1.05 you might need to tweak the lat and lon
+		     limits.
 
 
 ***CMIP5 ANALYSIS***
@@ -501,6 +502,33 @@ uv_vectorplot.py: Make a circumpolar Antarctic plot of speed overlaid with
 			  times as you like, and you only have to specify which
 			  parameters have changed since the last plot.
 
+massloss_map.py: Make a map of unexplained error in annually averaged simulated
+                 basal mass loss from each ice shelf that is over 5,000 km^2 in
+		 Rignot et al., 2013.
+		 To run: First make sure your logfile from massloss.py is up to
+		         date, as this script reads simulated basal mass loss
+			 timeseries from this logfile. Then open python or
+			 ipython and type "run massloss_map.py". The script will
+			 prompt you for paths to the ROMS grid file and the
+			 mass loss logfile, and whether to save the figure (and
+			 if so, what filename) or display it on the screen.
+
+nsidc_aice_monthly.py: Make a figure comparing sea ice concentration from NSIDC
+                       (1995 data) and CICE (latest year of spinup under
+		       repeated 1995 forcing), for the given month.
+		       Input: First edit the variables nsidc_head, nsidc_head_0,
+		              nsidc_head_1, and nsidc_tail (near the top of the
+			      file) to suit the paths to NSIDC output on your
+			      filesystem. Then open python or ipython and type
+			      "run nsidc_aice_monthly.py". The script will
+			      prompt you for the CICE output file, the month
+			      to plot, and whether to save the figure (and if
+			      so, what filename) or display it on the screen.
+			      As with the other interactive plotting scripts,
+			      the interface will repeat as many times as you
+			      like, and you only have to specify which
+			      parameters have changed since the last plot.
+
 
 ***ANIMATIONS***
 

massloss.py

@@ -141,15 +141,15 @@ def massloss (file_path, log_path):
         # Make sure y-limits won't cut off observed melt rate
         ymin = amin([ismr_low/factors[index], massloss_low, amin(massloss[index,:])])
         ymax = amax([ismr_high/factors[index], massloss_high, amax(massloss[index,:])])
-        if ymin < 0:
-            ymin = 1.05*ymin
-        else:
-            ymin = 0.95*ymin
-        if ymax < 0:
-            ymax = 0.95*ymax
-        else:
-            ymax = 1.05*ymax
-        ax1.set_ylim([ymin, ymax])
+        ticks = ax1.get_yticks()
+        min_tick = ticks[0]
+        max_tick = ticks[-1]
+        dtick = ticks[1]-ticks[0]
+        while min_tick >= ymin:
+            min_tick -= dtick
+        while max_tick <= ymax:
+            max_tick += dtick
+        ax1.set_ylim([min_tick, max_tick])
         # Title and ticks in blue for this side of the plot
         ax1.set_ylabel('Basal Mass Loss (Gt/y)', color='b')
         for t1 in ax1.get_yticklabels():

massloss_map.py

@@ -0,0 +1,175 @@
+from netCDF4 import Dataset
+from numpy import *
+from matplotlib.pyplot import *
+
+# Make a map of unexplained error in annually averaged simulated basal mass loss
+# from each ice shelf that is over 5,000 km^2 in Rignot et al., 2013.
+# Input:
+# grid_path = path to ROMS grid file
+# log_path = path to log file created by massloss.py
+# save = optional boolean to save the figure to a file, rather than displaying
+#        it on the screen
+# fig_name = if save=True, path to the desired filename for the figure
+def massloss_map (grid_path, log_path, save=False, fig_name=None):
+
+    # 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, -180, 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, 180]
+    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]
+
+    # Degrees to radians conversion factor
+    deg2rad = pi/180
+    # Northern boundary 60S 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
+
+    # Read log file
+    time = []
+    f = open(log_path, 'r')
+    # Skip the first line (header for time array)
+    f.readline()
+    for line in f:
+        try:
+            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
+    massloss_ts = empty([len(obs_massloss), len(time)])
+    index = 0
+    # Loop over ice shelves
+    while index < len(obs_massloss):
+        t = 0
+        for line in f:
+            try:
+                massloss_ts[index, t] = float(line)
+                t += 1
+            except(ValueError):
+                # Reached the header for the next variable
+                break
+        index +=1
+
+    # Find the time indices we care about: last year of simulation
+    time = array(time)
+    min_t = nonzero(time >= time[-1]-1)[0][0]
+    max_t = size(time)
+    # Find the average mass loss for each ice shelf over this last year
+    massloss = empty(len(obs_massloss))
+    for index in range(len(obs_massloss)):
+        massloss[index] = mean(massloss_ts[index, min_t:max_t])
+
+    # Read the grid
+    id = Dataset(grid_path, 'r')
+    lon = id.variables['lon_rho'][:,:]
+    lat = id.variables['lat_rho'][:,:]
+    mask_rho = id.variables['mask_rho'][:,:]
+    mask_zice = id.variables['mask_zice'][:,:]
+    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 a field of ice shelf mass loss unexplained error
+    error = zeros(shape(lon))
+    # Mask out land and open ocean points
+    error = ma.masked_where(mask_zice==0, error)
+
+    # Loop over ice shelves
+    for index in range(len(obs_massloss)):
+        # 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 error in mass loss
+        if massloss[index] < massloss_low:
+            # Simulated mass loss too low
+            error_tmp = massloss[index] - massloss_low
+        elif massloss[index] > massloss_high:
+            # Simulated mass loss too high
+            error_tmp = massloss[index] - massloss_high
+        else:
+            # Simulated mass loss within observational error estimates
+            error_tmp = 0
+        if error_tmp != 0:
+            # Modify error field for this region
+            if index == len(obs_massloss)-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)
+            error[region] = error_tmp
+
+    # 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)
+
+    # Determine bounds on colour scale
+    max_val = amax(abs(error))
+    lev = linspace(-max_val, max_val, num=40)
+
+    # Set up 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 error in mass loss
+    contourf(x, y, error, lev, cmap='RdBu_r')    
+    cbar = colorbar()
+    cbar.ax.tick_params(labelsize=20)
+    xlim([-nbdry, nbdry])
+    ylim([-nbdry, nbdry])
+    title('Bias in Ice Shelf Mass Loss (Gt/y)', fontsize=30)
+    axis('off')
+
+    # Finished
+    if save:
+        fig.savefig(fig_name)
+    else:
+        fig.show()
+
+
+# Command-line interface
+if __name__ == "__main__":
+
+    grid_path = raw_input("Path to grid file: ")
+    log_path = raw_input("Path to 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
+    # Make the plot
+    massloss_map(grid_path, log_path, save, fig_name)
+
+
+
+