Merge branch 'jhkennedy/stage' into 'master'

LEX-0.1-alpha

See merge request LIVVkit/lex-dev!2

clouds/clouds/cesm_cldhgh_ann.py

+import os
+import Ngl
+import numpy as np
+
+from netCDF4 import Dataset
+
+from livvkit.util import elements as el
+
+
+describe = """CESM_cldhgh_ANN plot."""
+
+
+def make_plot(config=None, out_path='.',
+              cloud_path='/lustre/atlas1/cli115/world-shared/4ue/obs_data/',
+              cism_path='/lustre/atlas1/cli115/world-shared/4ue/',
+              cesm_path='/lustre/atlas1/cli115/world-shared/4ue/b.e10.BG20TRCN.f09_g16.002/'):
+    # ---------------- Data source in TITAN ------------------------
+    f_cism = os.path.join(cism_path, 'Greenland_5km_v1.1_SacksRev_c110629.nc')
+    f_cesm_atm_climo_jja = os.path.join(cesm_path, 'postproc/atm/climos/b.e10.BG20TRCN.f09_g16.002_JJA_climo.nc')
+    f_isccp = os.path.join(cloud_path, 'ISCCP_ANN_climo.nc')
+    f_cloudsat = os.path.join(cloud_path, 'CLOUDSAT_ANN_climo.nc')
+    # --------------------------------------------------------------
+
+    img_list = []
+
+    # read f_cism, the elevation data
+    cism = Dataset(f_cism)
+    cism_usrf = cism.variables['usrf'][0, :, :]
+    cism_lat = cism.variables['lat'][0, :, :]
+    cism_lon = cism.variables['lon'][0, :, :]
+
+    # read f_cesm_atm_climo_jja and file2, CESM variable
+    cesm = Dataset(f_cesm_atm_climo_jja)
+    cesm_cldhgh = cesm.variables['CLDHGH'][0, :, :]
+    cesm_lat = cesm.variables['lat'][:]
+    cesm_lon = cesm.variables['lon'][:]
+
+    cesm_cldhgh = cesm_cldhgh * 100  # to get percent
+
+    # --- ISCCP data without remapping [64,128]
+    isccp = Dataset(f_isccp)
+    isccp_hgh = isccp.variables['CLDHGH'][0, :, :]
+    isccp_lat = isccp.variables['lat'][:]
+    isccp_lon = isccp.variables['lon'][:]
+
+    # --- CLOUDSAT data without remapping [94,192]
+    cloudsat = Dataset(f_cloudsat)
+    cloudsat_hgh = cloudsat.variables['CLDHGH'][:, :]
+    cloudsat_lat = cloudsat.variables['lat'][:]
+    cloudsat_lon = cloudsat.variables['lon'][:]
+
+    # print("contour plot of CLDHGH vs obs")
+    # maxmodel = np.max(cesm_cldhgh)
+    # maxisccp = np.max(isccp_hgh)
+    # maxcldsat = np.max(cloudsat_hgh)
+    # print("Max CLDHGH model: {}".format(maxmodel))
+    # print("Max CLDHGH f_isccp: {}".format(maxisccp))
+    # print("Max CLDHGH cldsat: {}".format(maxcldsat))
+
+    # ------- PLOT --------
+    # Open a workstation for drawing the plots
+    wkres = Ngl.Resources()
+    wkres.wkColorMap = "BlueWhiteOrangeRed"
+    # wkres.wkOrientation = "portrait"
+    wks_type = "png"
+    wks_img = str(os.path.join(out_path, 'CESM_cldhgh_ANN'))
+
+    wks = Ngl.open_wks(wks_type, wks_img, wkres)
+
+    # --- for the map -------
+    # Define plotting area, Greenland
+    mres = Ngl.Resources()
+    mres.nglDraw = False  # Don't draw individual plots
+    mres.nglFrame = False  # Don't advance frame.
+
+    mres.pmTickMarkDisplayMode = "Never"  # Turn off map tickmarks.
+    mres.mpGridAndLimbOn = False  # Turn off grid and limb lines.
+    mres.mpProjection = "Aitoff"
+    mres.mpLimitMode = "LatLon"  # limit map via lat/lon, to zoom in
+    mres.mpCenterLatF = 70.  # map area
+    mres.mpCenterLonF = -44.
+    mres.mpMinLatF = 57.
+    mres.mpMaxLatF = 85.
+    mres.mpMinLonF = -55.
+    mres.mpMaxLonF = -30.
+    mres.mpOutlineOn = False
+    mres.mpFillOn = False
+    mres.mpPerimOn = True  # add box around map
+
+    # --- for the CESM contour -------
+    res1 = Ngl.Resources()
+    res1.nglDraw = False  # Don't draw individual plots
+    res1.nglFrame = False  # Don't advance frame.
+    res1.cnLineLabelsOn = False
+    res1.cnFillOn = True
+    res1.cnLinesOn = False
+    res1.cnLineLabelsOn = False
+    res1.cnFillMode = "RasterFill"
+    res1.trGridType = "TriangularMesh"
+    res1.cnLevelSelectionMode = "ExplicitLevels"
+    res1.cnLevels = np.arange(0, 100, 10)
+    res1.lbLabelBarOn = True  # Turn on labelbar.
+    res1.lbLabelFontHeightF = 0.04
+
+    res1.sfXArray = cesm_lon
+    res1.sfYArray = cesm_lat
+
+    # --- for the data contour -------
+    res2 = Ngl.Resources()
+    res2.nglDraw = False  # Don't draw individual plots
+    res2.nglFrame = False  # Don't advance frame.
+    res2.cnLineLabelsOn = False
+    res2.cnFillOn = True
+    res2.cnLinesOn = False
+    res2.cnLineLabelsOn = False
+    res2.cnFillMode = "RasterFill"
+    res2.trGridType = "TriangularMesh"
+    res2.cnLevelSelectionMode = "ExplicitLevels"
+    res2.cnLevels = np.arange(0, 100, 10)
+    res2.lbLabelBarOn = True  # Turn on labelbar.
+    res2.lbOrientation = "Vertical"  # Verticle labelbar
+    res2.lbLabelFontHeightF = 0.04  # Make fonts smaller.
+
+    res2.sfXArray = isccp_lon
+    res2.sfYArray = isccp_lat
+
+    # --- for the data contour -------
+    res3 = Ngl.Resources()
+    res3.nglDraw = False  # Don't draw individual plots
+    res3.nglFrame = False  # Don't advance frame.
+    res3.cnLineLabelsOn = False
+    res3.cnFillOn = True
+    res3.cnLinesOn = False
+    res3.cnLineLabelsOn = False
+    res3.cnFillMode = "RasterFill"
+    res3.trGridType = "TriangularMesh"
+    res3.cnLevelSelectionMode = "ExplicitLevels"
+    res3.cnLevels = np.arange(0, 100, 10)
+    res3.lbLabelBarOn = True  # Turn on labelbar.
+    res3.lbOrientation = "Vertical"  # Verticle labelbar
+    res3.lbLabelFontHeightF = 0.04  # Make fonts smaller.
+
+    res3.sfXArray = cloudsat_lon
+    res3.sfYArray = cloudsat_lat
+
+    # ---- for the elevation -------
+    sres = Ngl.Resources()
+    sres.nglDraw = False  # Don't draw individual plots
+    sres.nglFrame = False  # Don't advance frame.
+    sres.cnFillOn = False
+    sres.cnLinesOn = True
+    sres.cnLineLabelsOn = False
+    sres.cnLevelSelectionMode = "ExplicitLevels"
+    sres.cnLevels = [0, 1000, 2000, 3000]
+    sres.sfXArray = cism_lon
+    sres.sfYArray = cism_lat
+
+    # ---- Overlay plots, each one has its own ID
+    # overlay ice on base map, and then overlay elevation on ice
+    usrf_plot1 = Ngl.contour(wks, cism_usrf, sres)
+    usrf_plot2 = Ngl.contour(wks, cism_usrf, sres)
+    usrf_plot3 = Ngl.contour(wks, cism_usrf, sres)
+
+    cesm_plot = Ngl.contour(wks, cesm_cldhgh, res1)
+    isccp_plot = Ngl.contour(wks, isccp_hgh, res2)
+    cldsat_plot = Ngl.contour(wks, cloudsat_hgh, res3)
+
+    # Create multiple figures and draw, which now contains the elevation and temperature
+    # "[1,3]" indicates 1 row, 3 columns.
+    map_title = ["CESM", "ISCCP", "CLOUDSAT"]
+
+    nmap = 3
+    plot = []
+    for i in range(nmap):
+        mres.tiMainString = map_title[i]
+        plot.append(Ngl.map(wks, mres))
+
+    # Overlay everything on the map plot.
+    Ngl.overlay(plot[0], cesm_plot)
+    Ngl.overlay(plot[0], usrf_plot1)
+    Ngl.overlay(plot[1], isccp_plot)
+    Ngl.overlay(plot[1], usrf_plot2)
+    Ngl.overlay(plot[2], cldsat_plot)
+    Ngl.overlay(plot[2], usrf_plot3)
+
+    Ngl.panel(wks, plot, [1, 3])
+
+    img_link = os.path.join(os.path.basename(out_path),
+                            os.path.basename(wks_img + '.' + wks_type))
+    img_elem = el.image('CESM_cldhgh_ANN',
+                        ' '.join(describe.split()),
+                        img_link)
+    if config:
+        img_elem['Height'] = config['image_height']
+    img_list.append(img_elem)
+
+    return img_list
+
+
+if __name__ == '__main__':
+    make_plot()

clouds/clouds/cesm_cldlow_ann.py

+import os
+import Ngl
+import numpy as np
+
+from netCDF4 import Dataset
+
+from livvkit.util import elements as el
+
+
+describe = """CESM_cldlow_ANN plot."""
+
+
+def make_plot(config=None, out_path='.',
+              cloud_path='/lustre/atlas1/cli115/world-shared/4ue/obs_data/',
+              cism_path='/lustre/atlas1/cli115/world-shared/4ue/',
+              cesm_path='/lustre/atlas1/cli115/world-shared/4ue/b.e10.BG20TRCN.f09_g16.002/'):
+    # ---------------- Data source in TITAN ------------------------
+    f_cism = os.path.join(cism_path, 'Greenland_5km_v1.1_SacksRev_c110629.nc')
+    f_cesm_atm_climo_jja = os.path.join(cesm_path, 'postproc/atm/climos/b.e10.BG20TRCN.f09_g16.002_JJA_climo.nc')
+    f_isccp = os.path.join(cloud_path, 'ISCCP_ANN_climo.nc')
+    f_cloudsat = os.path.join(cloud_path, 'CLOUDSAT_ANN_climo.nc')
+    # --------------------------------------------------------------
+
+    img_list = []
+
+    # read f_cism, the elevation data
+    cism = Dataset(f_cism)
+    cism_usrf = cism.variables['usrf'][0, :, :]
+    cism_lat = cism.variables['lat'][0, :, :]
+    cism_lon = cism.variables['lon'][0, :, :]
+
+    # read f_cesm_atm_climo_jja, CESM variable
+    cesm = Dataset(f_cesm_atm_climo_jja)
+    cesm_low = cesm.variables['CLDLOW'][0, :, :]
+    cesm_lat = cesm.variables['lat'][:]
+    cesm_lon = cesm.variables['lon'][:]
+
+    cesm_low = cesm_low * 100  # to get percent
+
+    # --- ISCCP data without remapping [64,128]
+    isccp = Dataset(f_isccp)
+    isccp_low = isccp.variables['CLDLOW'][0, :, :]
+    isccp_lat = isccp.variables['lat'][:]
+    isccp_lon = isccp.variables['lon'][:]
+
+    # --- CLOUDSAT data without remapping [94,192]
+    cloudsat = Dataset(f_cloudsat)
+    cloudsat_low = cloudsat.variables['CLDLOW'][:, :]
+    cloudsat_lat = cloudsat.variables['lat'][:]
+    cloudsat_lon = cloudsat.variables['lon'][:]
+
+    # print("contour plot of CLDLOW vs obs")
+    # maxmodel = np.max(cesm_low)
+    # maxisccp = np.max(isccp_low)
+    # maxcldsat = np.max(cloudsat_low)
+    # print("Max CLDLOW model:  {}".format(maxmodel))
+    # print("Max CLDLOW isccp:  {}".format(maxisccp))
+    # print("Max CLDLOW cldsat: {}".format(maxcldsat))
+
+    # ------- PLOT --------
+    #  Open a workstation for drawing the plots
+    wkres = Ngl.Resources()
+    # wkres.wkColorMap = "WhiteBlueGreenYellowRed"
+    wkres.wkColorMap = "BlueWhiteOrangeRed"
+    # wkres.wkOrientation = "portrait"  # "portrait" or "landscape"
+    wks_type = "png"
+    wks_img = str(os.path.join(out_path, "CESM_cldlow_ANN"))
+    wks = Ngl.open_wks(wks_type, wks_img, wkres)
+
+    # --- for the map -------
+    # Define plotting area, Greenland
+    mres = Ngl.Resources()
+    mres.nglDraw = False  # Don't draw individual plots
+    mres.nglFrame = False  # Don't advance frame.
+
+    mres.pmTickMarkDisplayMode = "Never"  # Turn off map tickmarks.
+    mres.mpGridAndLimbOn = False  # Turn off grid and limb lines.
+    mres.mpProjection = "Aitoff"
+    mres.mpLimitMode = "LatLon"  # limit map via lat/lon, to zoom in
+    mres.mpCenterLatF = 70.  # map area
+    mres.mpCenterLonF = -44.
+    mres.mpMinLatF = 57.
+    mres.mpMaxLatF = 85.
+    mres.mpMinLonF = -55.
+    mres.mpMaxLonF = -30.
+    mres.mpOutlineOn = False
+    mres.mpFillOn = False
+    mres.mpPerimOn = True  # add box around map
+
+    # --- for the CESM contour -------
+    res1 = Ngl.Resources()
+    res1.nglDraw = False  # Don't draw individual plots
+    res1.nglFrame = False  # Don't advance frame.
+    res1.cnLineLabelsOn = False
+    res1.cnFillOn = True
+    res1.cnLinesOn = False
+    res1.cnLineLabelsOn = False
+    res1.cnFillMode = "RasterFill"
+    res1.trGridType = "TriangularMesh"
+    res1.cnLevelSelectionMode = "ExplicitLevels"
+    res1.cnLevels = np.arange(0, 100, 10)
+    res1.lbLabelBarOn = True  # Turn on labelbar.
+    res1.lbLabelFontHeightF = 0.04
+
+    res1.sfXArray = cesm_lon
+    res1.sfYArray = cesm_lat
+
+    # --- for the data contour -------
+    res2 = Ngl.Resources()
+    res2.nglDraw = False  # Don't draw individual plots
+    res2.nglFrame = False  # Don't advance frame.
+    res2.cnLineLabelsOn = False
+    res2.cnFillOn = True
+    res2.cnLinesOn = False
+    res2.cnLineLabelsOn = False
+    res2.cnFillMode = "RasterFill"
+    res2.trGridType = "TriangularMesh"
+    res2.cnLevelSelectionMode = "ExplicitLevels"
+    res2.cnLevels = np.arange(0, 100, 10)
+    res2.lbLabelBarOn = True  # Turn on labelbar.
+    res2.lbOrientation = "Vertical"  # Verticle labelbar
+    res2.lbLabelFontHeightF = 0.04  # Make fonts smaller.
+
+    res2.sfXArray = isccp_lon
+    res2.sfYArray = isccp_lat
+
+    # --- for the data contour -------
+    res3 = Ngl.Resources()
+    res3.nglDraw = False  # Don't draw individual plots
+    res3.nglFrame = False  # Don't advance frame.
+    res3.cnLineLabelsOn = False
+    res3.cnFillOn = True
+    res3.cnLinesOn = False
+    res3.cnLineLabelsOn = False
+    res3.cnFillMode = "RasterFill"
+    res3.trGridType = "TriangularMesh"
+    res3.cnLevelSelectionMode = "ExplicitLevels"
+    res3.cnLevels = np.arange(0, 100, 10)
+    res3.lbLabelBarOn = True  # Turn on labelbar.
+    res3.lbOrientation = "Vertical"  # Verticle labelbar
+    res3.lbLabelFontHeightF = 0.04  # Make fonts smaller.
+
+    res3.sfXArray = cloudsat_lon
+    res3.sfYArray = cloudsat_lat
+
+    # ---- for the elevation -------
+    sres = Ngl.Resources()
+    sres.nglDraw = False  # Don't draw individual plots
+    sres.nglFrame = False  # Don't advance frame.
+    sres.cnFillOn = False
+    sres.cnLinesOn = True
+    sres.cnLineLabelsOn = False
+    sres.cnLevelSelectionMode = "ExplicitLevels"
+    sres.cnLevels = [0, 1000, 2000, 3000]
+    sres.sfXArray = cism_lon
+    sres.sfYArray = cism_lat
+
+    # ---- Overlay plots, each one has its own ID
+    # overlay ice on base map, and then overlay elevation on ice
+    usrf_plot1 = Ngl.contour(wks, cism_usrf, sres)
+    usrf_plot2 = Ngl.contour(wks, cism_usrf, sres)
+    usrf_plot3 = Ngl.contour(wks, cism_usrf, sres)
+
+    cesm_plot = Ngl.contour(wks, cesm_low, res1)
+    isccp_plot = Ngl.contour(wks, isccp_low, res2)
+    cldsat_plot = Ngl.contour(wks, cloudsat_low, res3)
+
+    # Creat multiple figures and draw, which now contains the elevation and temperature
+    # "[1,3]" indicates 1 row, 3 columns.
+    map_title = ["CESM", "ISCCP", "CLOUDSAT"]
+
+    nmap = 3
+    plot = []
+    for i in range(nmap):
+        mres.tiMainString = map_title[i]
+        plot.append(Ngl.map(wks, mres))
+
+    # Overlay everything on the map plot.
+    Ngl.overlay(plot[0], cesm_plot)
+    Ngl.overlay(plot[0], usrf_plot1)
+    Ngl.overlay(plot[1], isccp_plot)
+    Ngl.overlay(plot[1], usrf_plot2)
+    Ngl.overlay(plot[2], cldsat_plot)
+    Ngl.overlay(plot[2], usrf_plot3)
+
+    Ngl.panel(wks, plot, [1, 3])
+
+    img_link = os.path.join(os.path.basename(out_path),
+                            os.path.basename(wks_img + '.' + wks_type))
+    img_elem = el.image('CESM_cldlow_ANN',
+                        ' '.join(describe.split()),
+                        img_link)
+    if config:
+        img_elem['Height'] = config['image_height']
+    img_list.append(img_elem)
+
+    return img_list
+
+
+if __name__ == '__main__':
+    make_plot()

clouds/clouds/cesm_cldtot_ann.py

+import os
+import Ngl
+import numpy as np
+
+from netCDF4 import Dataset
+
+from livvkit.util import elements as el
+
+
+describe = """CESM_cldtot_ANN plot."""
+
+
+def make_plot(config=None, out_path='.',
+              cloud_path='/lustre/atlas1/cli115/world-shared/4ue/obs_data/',
+              cism_path='/lustre/atlas1/cli115/world-shared/4ue/',
+              cesm_path='/lustre/atlas1/cli115/world-shared/4ue/b.e10.BG20TRCN.f09_g16.002/'):
+    # ---------------- Data source in TITAN ------------------------
+    f_cism = os.path.join(cism_path, 'Greenland_5km_v1.1_SacksRev_c110629.nc')
+    f_cesm_atm_climo_jja = os.path.join(cesm_path, 'postproc/atm/climos/b.e10.BG20TRCN.f09_g16.002_JJA_climo.nc')
+    f_isccp = os.path.join(cloud_path, 'ISCCP_ANN_climo.nc')
+    f_cloudsat = os.path.join(cloud_path, 'CLOUDSAT_ANN_climo.nc')
+    # --------------------------------------------------------------
+
+    img_list = []
+
+    # read f_cism, the elevation data
+    cism = Dataset(f_cism)
+    cism_usrf = cism.variables['usrf'][0, :, :]
+    cism_lat = cism.variables['lat'][0, :, :]
+    cism_lon = cism.variables['lon'][0, :, :]
+    
+    # read f_cesm_atm_climo_jja and file2, CESM variable
+    cesm = Dataset(f_cesm_atm_climo_jja)
+    cesm_tot = cesm.variables['CLDTOT'][0, :, :]
+    cesm_lat = cesm.variables['lat'][:]
+    cesm_lon = cesm.variables['lon'][:]
+    
+    cesm_tot = cesm_tot * 100  # to get percent
+    
+    # --- ISCCP data without remapping [64,128]
+    isccp = Dataset(f_isccp)
+    isccp_tot = isccp.variables['CLDTOT'][0, :, :]
+    isccp_lat = isccp.variables['lat'][:]
+    isccp_lon = isccp.variables['lon'][:]
+    
+    # --- CLOUDSAT data without remapping [94,192]
+    cloudsat = Dataset(f_cloudsat)
+    cloudsat_tot = cloudsat.variables['CLDTOT'][:, :]
+    cloudsat_lat = cloudsat.variables['lat'][:]
+    cloudsat_lon = cloudsat.variables['lon'][:]
+    
+    # print("contour plot of CLDTOT vs obs")
+    # maxmodel = np.max(cesm_tot)
+    # maxisccp = np.max(isccp_tot)
+    # maxcldsat = np.max(cloudsat_tot)
+    # print("Max CLDTOT model:  {}".format(maxmodel))
+    # print("Max CLDTOT isccp:  {}".format(maxisccp))
+    # print("Max CLDTOT cldsat: {}".format(maxcldsat))
+    
+    # ------- PLOT --------
+    #  Open a workstation for drawing the plots
+    wkres = Ngl.Resources()
+    # wkres.wkColorMap = "WhiteBlueGreenYellowRed"
+    wkres.wkColorMap = "BlueWhiteOrangeRed"
+    # wkres.wkOrientation = "portrait"  # "portrait" or "landscape"
+    wks_type = "png"
+    wks_img = str(os.path.join(out_path, "CESM_cldtot_ANN"))
+    wks = Ngl.open_wks(wks_type, wks_img, wkres)
+    
+    # --- for the map -------
+    # Define plotting area, Greenland
+    mres = Ngl.Resources()
+    mres.nglDraw = False  # Don't draw individual plots
+    mres.nglFrame = False  # Don't advance frame.
+    
+    mres.pmTickMarkDisplayMode = "Never"  # Turn off map tickmarks.
+    mres.mpGridAndLimbOn = False  # Turn off grid and limb lines.
+    mres.mpProjection = "Aitoff"
+    mres.mpLimitMode = "LatLon"  # limit map via lat/lon, to zoom in
+    mres.mpCenterLatF = 70.  # map area
+    mres.mpCenterLonF = -44.
+    mres.mpMinLatF = 57.
+    mres.mpMaxLatF = 85.
+    mres.mpMinLonF = -55.
+    mres.mpMaxLonF = -30.
+    mres.mpOutlineOn = False
+    mres.mpFillOn = False
+    mres.mpPerimOn = True  # add box around map
+    
+    # --- for the CESM contour -------
+    res1 = Ngl.Resources()
+    res1.nglDraw = False  # Don't draw individual plots
+    res1.nglFrame = False  # Don't advance frame.
+    res1.cnLineLabelsOn = False
+    res1.cnFillOn = True
+    res1.cnLinesOn = False
+    res1.cnLineLabelsOn = False
+    res1.cnFillMode = "RasterFill"
+    res1.trGridType = "TriangularMesh"
+    res1.cnLevelSelectionMode = "ExplicitLevels"
+    res1.cnLevels = np.arange(0, 100, 10)
+    res1.lbLabelBarOn = True  # Turn on labelbar.
+    res1.lbLabelFontHeightF = 0.04
+    
+    res1.sfXArray = cesm_lon
+    res1.sfYArray = cesm_lat
+    
+    # --- for the data contour -------
+    res2 = Ngl.Resources()
+    res2.nglDraw = False  # Don't draw individual plots
+    res2.nglFrame = False  # Don't advance frame.
+    res2.cnLineLabelsOn = False
+    res2.cnFillOn = True
+    res2.cnLinesOn = False
+    res2.cnLineLabelsOn = False
+    res2.cnFillMode = "RasterFill"
+    res2.trGridType = "TriangularMesh"
+    res2.cnLevelSelectionMode = "ExplicitLevels"
+    res2.cnLevels = np.arange(0, 100, 10)
+    res2.lbLabelBarOn = True  # Turn on labelbar.
+    res2.lbOrientation = "Vertical"  # Verticle labelbar
+    res2.lbLabelFontHeightF = 0.04  # Make fonts smaller.
+    
+    res2.sfXArray = isccp_lon
+    res2.sfYArray = isccp_lat
+    
+    # --- for the data contour -------
+    res3 = Ngl.Resources()
+    res3.nglDraw = False  # Don't draw individual plots
+    res3.nglFrame = False  # Don't advance frame.
+    res3.cnLineLabelsOn = False
+    res3.cnFillOn = True
+    res3.cnLinesOn = False
+    res3.cnLineLabelsOn = False
+    res3.cnFillMode = "RasterFill"
+    res3.trGridType = "TriangularMesh"
+    res3.cnLevelSelectionMode = "ExplicitLevels"
+    res3.cnLevels = np.arange(0, 100, 10)
+    res3.lbLabelBarOn = True  # Turn on labelbar.
+    res3.lbOrientation = "Vertical"  # Verticle labelbar
+    res3.lbLabelFontHeightF = 0.04  # Make fonts smaller.
+    
+    res3.sfXArray = cloudsat_lon
+    res3.sfYArray = cloudsat_lat
+    
+    # ---- for the elevation -------
+    sres = Ngl.Resources()
+    sres.nglDraw = False  # Don't draw individual plots
+    sres.nglFrame = False  # Don't advance frame.
+    sres.cnFillOn = False
+    sres.cnLinesOn = True
+    sres.cnLineLabelsOn = False
+    sres.cnLevelSelectionMode = "ExplicitLevels"
+    sres.cnLevels = [0, 1000, 2000, 3000]
+    sres.sfXArray = cism_lon
+    sres.sfYArray = cism_lat
+    
+    # ---- Overlay plots, each one has its own ID
+    # overlay ice on base map, and then overlay elevation on ice
+    usrf_plot1 = Ngl.contour(wks, cism_usrf, sres)
+    usrf_plot2 = Ngl.contour(wks, cism_usrf, sres)
+    usrf_plot3 = Ngl.contour(wks, cism_usrf, sres)
+    
+    cesm_plot = Ngl.contour(wks, cesm_tot, res1)
+    isccp_plot = Ngl.contour(wks, isccp_tot, res2)
+    cldsat_plot = Ngl.contour(wks, cloudsat_tot, res3)
+    
+    # Creat multiple figures and draw, which now contains the elevation and temperature
+    # "[1,3]" indicates 1 row, 3 columns.
+    map_title = ["CESM", "ISCCP", "CLOUDSAT"]
+    
+    nmap = 3
+    plot = []
+    for i in range(nmap):
+        mres.tiMainString = map_title[i]
+        plot.append(Ngl.map(wks, mres))
+    
+    # Overlay everything on the map plot.
+    Ngl.overlay(plot[0], cesm_plot)
+    Ngl.overlay(plot[0], usrf_plot1)
+    Ngl.overlay(plot[1], isccp_plot)
+    Ngl.overlay(plot[1], usrf_plot2)
+    Ngl.overlay(plot[2], cldsat_plot)
+    Ngl.overlay(plot[2], usrf_plot3)
+    
+    Ngl.panel(wks, plot, [1, 3])
+
+    img_link = os.path.join(os.path.basename(out_path),
+                            os.path.basename(wks_img + '.' + wks_type))