import os import Ngl import numpy as np from netCDF4 import Dataset from livvkit.util import elements as el describe = """ Climatological annual average of low cloud cover over Greenland for CESM (left), ISCCP (middle; Rossow and Schiffer, 1999), and CLOUDSAT (right; Kay and Gettelman, 2009). The black solid lines denote the 0, 1000, 2000, and 3000 meter Greenland ice sheet elevation contours as seen by CESM. """ title = "Low cloud cover over Greenland" def make_plot(config, out_path='.'): # ---------------- Data source in TITAN ------------------------ f_cism = os.path.join(config['cism_data'], 'Greenland_5km_v1.1_SacksRev_c110629.nc') f_cesm_atm_climo_jja = os.path.join(config['cesm_atm_climos'], 'b.e10.BG20TRCN.f09_g16.002_JJA_climo.nc') f_isccp = os.path.join(config['cloud_data'], 'ISCCP_ANN_climo.nc') f_cloudsat = os.path.join(config['cloud_data'], '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 #cmap = ["grey98","grey94","grey90","grey86","grey80","grey74","grey68","grey60","grey52","grey44","grey36"] wkres = Ngl.Resources() # 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.cnFillPalette = "percent_11lev" 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.cnFillPalette = "percent_11lev" 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.cnFillPalette = "percent_11lev" 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(title, ' '.join(describe.split()), img_link) img_elem['Height'] = config['image_height'] img_list.append(img_elem) return img_list