More intercomparison figures

mip_drift_slices.py

@@ -11,11 +11,12 @@
 sys.path.insert(0, '/short/y99/kaa561/fesomtools')
 from fesom_grid import *
 from fesom_sidegrid import *
+from triangle_area import *
 
 # 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.
+# are the last January of the simulation (monthly average) from MetROMS and
+# FESOM respectively.
 # Input:
 # roms_grid = path to ROMS grid file
 # roms_file = path to file containing Jan 2016 monthly average of temperature
@@ -45,6 +46,14 @@ def mip_drift_slices (roms_grid, roms_file, fesom_mesh_path, fesom_file):
     temp_max = 6
     salt_min = 33.9
     salt_max = 34.9
+    # Contours to overlay
+    temp_contour = 0.75
+    salt_contour = 34.5
+    # Parameters for FESOM regular grid interpolation (needed for contours)
+    num_lat = 1000
+    num_depth = 500
+    r = 6.371e6
+    deg2rad = pi/180.0
 
     # Get longitude for the title
     if lon0 < 0:
@@ -67,7 +76,7 @@ def mip_drift_slices (roms_grid, roms_file, fesom_mesh_path, fesom_file):
         id.close()
     else:
         print 'lon0 is only coded for 0E at this time'
-        return
+        #return
 
     print 'Processing ROMS'
     # Read grid variables we need
@@ -118,7 +127,75 @@ def mip_drift_slices (roms_grid, roms_file, fesom_mesh_path, fesom_file):
     # Repeat for the other variables
     fesom_salt = []
     for selm in selements_salt:
-        fesom_salt.append(selm.var)    
+        fesom_salt.append(selm.var)
+    # Interpolate to regular grid so we can overlay contours
+    lat_reg = linspace(lat_min, lat_max, num_lat)
+    depth_reg = linspace(-depth_max, -depth_min, num_depth)
+    temp_reg = zeros([num_depth, num_lat])
+    salt_reg = zeros([num_depth, num_lat])
+    temp_reg[:,:] = NaN
+    salt_reg[:,:] = NaN
+    # For each element, check if a point on the regular grid lies
+    # within. If so, do barycentric interpolation to that point, at each
+    # depth on the regular grid.
+    for elm in elements:
+        # Check if this element crosses lon0
+        if amin(elm.lon) < lon0 and amax(elm.lon) > lon0:
+            # Check if we are within the latitude bounds
+            if amax(elm.lat) > lat_min and amin(elm.lat) < lat_max:
+                # Find largest regular latitude value south of Element
+                tmp = nonzero(lat_reg > amin(elm.lat))[0]
+                if len(tmp) == 0:
+                    # Element crosses the southern boundary
+                    jS = 0
+                else:
+                    jS = tmp[0] - 1
+                # Find smallest regular latitude north of Element
+                tmp = nonzero(lat_reg > amax(elm.lat))[0]
+                if len(tmp) == 0:
+                    # Element crosses the northern boundary
+                    jN = num_lat
+                else:
+                    jN = tmp[0]
+                for j in range(jS+1,jN):
+                    # There is a chance that the regular gridpoint at j
+                    # lies within this element
+                    lat0 = lat_reg[j]
+                    if in_triangle(elm, lon0, lat0):
+                        # Yes it does
+                        # Get area of entire triangle
+                        area = triangle_area(elm.lon, elm.lat)
+                        # Get area of each sub-triangle formed by (lon0, lat0)
+                        area0 = triangle_area([lon0, elm.lon[1], elm.lon[2]], [lat0, elm.lat[1], elm.lat[2]])
+                        area1 = triangle_area([lon0, elm.lon[0], elm.lon[2]], [lat0, elm.lat[0], elm.lat[2]])
+                        area2 = triangle_area([lon0, elm.lon[0], elm.lon[1]], [lat0, elm.lat[0], elm.lat[1]])
+                        # Find fractional area of each
+                        cff = [area0/area, area1/area, area2/area]
+                        # Interpolate each depth value
+                        for k in range(num_depth):
+                            # Linear interpolation in the vertical for the
+                            # value at each corner of the triangle
+                            node_vals_temp = []
+                            node_vals_salt = []
+                            for n in range(3):
+                                id1, id2, coeff1, coeff2 = elm.nodes[n].find_depth(depth_reg[k])
+                                if any(isnan(array([id1, id2, coeff1, coeff2]))):
+                                    # No ocean data here (seafloor or ice shelf)
+                                    node_vals_temp.append(NaN)
+                                    node_vals_salt.append(NaN)
+                                else:
+                                    node_vals_temp.append(coeff1*fesom_temp_nodes[id1] + coeff2*fesom_temp_nodes[id2])
+                                    node_vals_salt.append(coeff1*fesom_salt_nodes[id1] + coeff2*fesom_salt_nodes[id2])
+                            if any(isnan(node_vals_temp)):
+                                pass
+                            else:
+                                # Barycentric interpolation for the value at
+                                # lon0, lat0
+                                temp_reg[k,j] = sum(array(cff)*array(node_vals_temp))
+                                salt_reg[k,j] = sum(array(cff)*array(node_vals_salt))
+    temp_reg = ma.masked_where(isnan(temp_reg), temp_reg)
+    salt_reg = ma.masked_where(isnan(salt_reg), salt_reg)
+    depth_reg = -1*depth_reg
 
     # Set up axis labels the way we want them
     lat_ticks = arange(lat_min+3, lat_max+10, 10)
@@ -138,6 +215,8 @@ def mip_drift_slices (roms_grid, roms_file, fesom_mesh_path, fesom_file):
     # Temperature
     ax = subplot(gs1[0,0])
     pcolor(ecco_lat, ecco_depth, ecco_temp, vmin=temp_min, vmax=temp_max, cmap='jet')
+    # Overlay contour
+    contour(ecco_lat, ecco_depth, ecco_temp, levels=[temp_contour], color='black')
     title(r'Temperature ($^{\circ}$C)', fontsize=24)
     ylabel('Depth (m)', fontsize=18)
     xlim([lat_min, lat_max])
@@ -150,6 +229,7 @@ def mip_drift_slices (roms_grid, roms_file, fesom_mesh_path, fesom_file):
     # Salinity
     ax = subplot(gs1[0,1])
     pcolor(ecco_lat, ecco_depth, ecco_salt, vmin=salt_min, vmax=salt_max, cmap='jet')
+    contour(ecco_lat, ecco_depth, ecco_salt, levels=[salt_contour], color='black')
     title('Salinity (psu)', fontsize=24)
     xlim([lat_min, lat_max])
     ylim([depth_min, depth_max])
@@ -163,6 +243,7 @@ def mip_drift_slices (roms_grid, roms_file, fesom_mesh_path, fesom_file):
     # Temperature
     ax = subplot(gs2[0,0])
     pcolor(roms_lat, roms_z, roms_temp, vmin=temp_min, vmax=temp_max, cmap='jet')
+    contour(roms_lat, roms_z, roms_temp, levels=[temp_contour], color='black')
     ylabel('Depth (m)', fontsize=18)
     xlim([lat_min, lat_max])
     ylim([depth_min, depth_max])
@@ -174,6 +255,7 @@ def mip_drift_slices (roms_grid, roms_file, fesom_mesh_path, fesom_file):
     # Salinity
     ax = subplot(gs2[0,1])
     pcolor(roms_lat, roms_z, roms_salt, vmin=salt_min, vmax=salt_max, cmap='jet')
+    contour(roms_lat, roms_z, roms_salt, levels=[salt_contour], color='black')
     xlim([lat_min, lat_max])
     ylim([depth_min, depth_max])
     ax.set_xticks(lat_ticks)
@@ -190,6 +272,8 @@ def mip_drift_slices (roms_grid, roms_file, fesom_mesh_path, fesom_file):
     img.set_edgecolor('face')
     img.set_clim(vmin=temp_min, vmax=temp_max)
     ax.add_collection(img)
+    # Overlay contour on regular grid
+    contour(lat_reg, depth_reg, temp_reg, levels=[temp_contour], color='black')
     ylabel('Depth (m)', fontsize=18)
     xlim([lat_min, lat_max])
     ylim([depth_min, depth_max])
@@ -209,6 +293,7 @@ def mip_drift_slices (roms_grid, roms_file, fesom_mesh_path, fesom_file):
     img.set_edgecolor('face')
     img.set_clim(vmin=salt_min, vmax=salt_max)
     ax.add_collection(img)
+    contour(lat_reg, depth_reg, salt_reg, levels=[salt_contour], color='black')
     xlim([lat_min, lat_max])
     ylim([depth_min, depth_max])
     ax.set_xticks(lat_ticks)
@@ -223,13 +308,22 @@ def mip_drift_slices (roms_grid, roms_file, fesom_mesh_path, fesom_file):
     fig.savefig('ts_drift.png')
 
 
+def in_triangle (elm, lon0, lat0):
+
+   alpha = ((elm.lat[1] - elm.lat[2])*(lon0 - elm.lon[2]) + (elm.lon[2] - elm.lon[1])*(lat0 - elm.lat[2]))/((elm.lat[1] - elm.lat[2])*(elm.lon[0] - elm.lon[2]) + (elm.lon[2] - elm.lon[1])*(elm.lat[0] - elm.lat[2]))
+   beta = ((elm.lat[2] - elm.lat[0])*(lon0 - elm.lon[2]) + (elm.lon[0] - elm.lon[2])*(lat0 - elm.lat[2]))/((elm.lat[1] - elm.lat[2])*(elm.lon[0] - elm.lon[2]) + (elm.lon[2] - elm.lon[1])*(elm.lat[0] - elm.lat[2]))
+   gamma = 1 - alpha - beta
+
+   return alpha >= 0 and beta >= 0 and gamma >= 0
+
+
 # Command-line interface
 if __name__ == "__main__":
 
     roms_grid = raw_input("Path to ROMS grid file: ")
-    roms_file = raw_input("Path to last ocean_avg file output by ROMS: ")
+    roms_file = raw_input("Path to ROMS file containing monthly averaged temperature and salinity for January 2016: ")
     fesom_mesh_path = raw_input("Path to FESOM mesh directory: ")
-    fesom_file = raw_input("Path to last oce.mean file output by FESOM: ")
+    fesom_file = raw_input("Path to FESOM file containing monthly averaged temperature and salinity for January 2016: ")
     mip_drift_slices(roms_grid, roms_file, fesom_mesh_path, fesom_file)
 
 

mip_scatterplot.py

@@ -0,0 +1,189 @@
+from numpy import *
+from matplotlib.pyplot import *
+
+def mip_scatterplot (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 = ['Larsen D', 'Larsen C', 'Wilkins & George VI & Stange', 'Filchner-Ronne', 'Abbot', 'Pine Island Glacier', 'Thwaites', 'Dotson', 'Getz', 'Nickerson', 'Sulzberger', 'Mertz', 'Totten & Moscow University', 'Shackleton', 'West', 'Amery', 'Prince Harald', 'Baudouin & Borchgrevink', 'Lazarev', 'Nivl', 'Fimbul & Jelbart & Ekstrom', 'Brunt & Riiser-Larsen', 'Ross']
+    # 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)
+    # Order of indices for the ice shelves to be plotted on the x-axis (0-based)
+    order = [3, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 22, 10, 9, 8, 7, 6, 5, 4, 2, 1, 0]
+
+    # 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)])
+    # Skip total mass loss
+    for line in f:
+        try:
+            tmp = float(line)
+        except(ValueError):
+            # Reaced the header for the first ice shelf
+            break
+    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
+    f = open(fesom_logfile_lr, 'r')
+    # Skip the first line (header)
+    f.readline()
+    # Skip total mass loss
+    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_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_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
+    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
+    f.close()
+    fesom_massloss_hr = mean(fesom_massloss_ts_hr[:,peryear*(calc_start-year_start):peryear*(calc_end+1-year_start)], axis=1)
+
+    # Figure out error values, in correct order for plotting
+    roms_error = []
+    fesom_error_lr = []
+    fesom_error_hr = []
+    labels = []
+    for index in range(num_shelves):
+        obs_min = obs_massloss[order[index]] - obs_massloss_error[order[index]]
+        obs_max = obs_massloss[order[index]] + obs_massloss_error[order[index]]
+        if roms_massloss[order[index]] < obs_min:
+            roms_error.append(roms_massloss[order[index]] - obs_min)
+        elif roms_massloss[order[index]] > obs_max:
+            roms_error.append(roms_massloss[order[index]] - obs_max)
+        else:
+            roms_error.append(0)
+        if fesom_massloss_lr[order[index]] < obs_min:
+            fesom_error_lr.append(fesom_massloss_lr[order[index]] - obs_min)
+        elif fesom_massloss_lr[order[index]] > obs_max:
+            fesom_error_lr.append(fesom_massloss_lr[order[index]] - obs_max)
+        else:
+            fesom_error_lr.append(0)
+        if fesom_massloss_hr[order[index]] < obs_min:
+            fesom_error_hr.append(fesom_massloss_hr[order[index]] - obs_min)
+        elif fesom_massloss_hr[order[index]] > obs_max:
+            fesom_error_hr.append(fesom_massloss_hr[order[index]] - obs_max)
+        else:
+            fesom_error_hr.append(0)
+        labels.append(names[order[index]])
+
+    # Plot
+    fig = figure(figsize=(10,7))
+    gs = GridSpec(1,1)
+    gs.update(left=0.1, right=0.9, bottom=0.4, top=0.9)
+    ax = subplot(gs[0,0])
+    # Alternate background between white and light blue to split up regions
+    axvspan(0.5, 6.5, facecolor='b', alpha=0.1)
+    axvspan(7.5, 11.5, facecolor='b', alpha=0.1)
+    axvspan(14.5, 18.5, facecolor='b', alpha=0.1)
+    axvspan(20.5, 23, facecolor='b', alpha=0.1)
+    # Region labels
+    text(0, 30, 'FR', fontsize=14, ha='center')
+    text(3.5, 30, 'EWed', fontsize=14, ha='center')
+    text(7, -20, 'Am', fontsize=14, ha='center')
+    text(9.5, 30, 'Aus', fontsize=14, ha='center')
+    text(13, 30, 'RS', fontsize=14, ha='center')
+    text(16.5, 30, 'AS', fontsize=14, ha='center')
+    text(19.5, 30, 'BS', fontsize=14, ha='center')
+    text(21.5, 30, 'Lr', fontsize=14, ha='center')
+    plot(range(-1,num_shelves+1), zeros(num_shelves+2), color='black', linewidth=3)
+    plot(range(num_shelves), roms_error, 'o', color=(0.08, 0.4, 0.79), ms=10, label='MetROMS')
+    plot(range(num_shelves), fesom_error_lr, 'o', color=(0.73, 0.06, 0.69), ms=10, label='FESOM (low-res)')
+    plot(range(num_shelves), fesom_error_hr, 'o', color=(0.06, 0.73, 0.1), ms=10, label='FESOM (high-res)')
+    grid(True)    
+    xlim([-0.5, num_shelves-0.5])
+    xticks(range(num_shelves), labels, rotation=90)
+    ylabel('Gt/y', fontsize=14)
+    title('Bias in Ice Shelf Basal Mass Loss', fontsize=18)
+    # Make legend
+    legend(numpoints=1,loc='lower left')
+    setp(gca().get_legend().get_texts(), fontsize='13')
+    fig.show()
+
+
+# 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_scatterplot(roms_logfile, fesom_logfile_lr, fesom_logfile_hr)
+