# coding=utf-8 import os import Ngl import numpy as np import numpy.ma as ma from netCDF4 import Dataset from livvkit.util import elements as el describe = """ Average of the annual average latent heat flux (W m^-2) over Greenland for every summer (June-July-August; JJA) from 1980--1999 for CESM (left) and RACMO 2.3 (middle), and the difference between them (right; CESM - RACMO 2.3). The black solid lines denote the 0, 1000, 2000, and 3000 meter Greenland ice sheet elevation contours as seen by the models (CESM's contours shown in the difference plot). """ title = "Summer latent heat flux" 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_lnd_climo_jja = os.path.join(config['cesm_lnd_climos'], 'b.e10.BG20TRCN.f09_g16.002_JJA_196006_200508_climo.nc') f_racmo_latf_jja = os.path.join(config['racmo_data'], 'racmo23_GRN_monthly.latf.1980-1999.JJA.nc') f_racmo_latf_jja_remapped = os.path.join(config['racmo_data'], 'racmo23_GRN_monthly.latf.1980-1999.remap2cesm.JJA.nc') f_racmo_mask = os.path.join(config['racmo_data'], 'RACMO23_masks_ZGRN11.nc') # -------------------------------------------------------------- img_list = [] # f_cism get the following, matrixsize[301,561] # usrf(0,:,:), lat(0,:,:), lon(0,:,:) # Goal: plot 3 figures, radiation overlaying on elevation usrf # read f_cism, elevation 5km CESM # adding the [:] transforms the variable to type numpy.ndarray ncid0 = Dataset(f_cism) usrf = ncid0.variables['usrf'][0, :, :] lat = ncid0.variables['lat'][0, :, :] lon = ncid0.variables['lon'][0, :, :] # read f_cesm_lnd_climo_jja, CESM variables ncid1 = Dataset(f_cesm_lnd_climo_jja) qsoil = ncid1.variables['QSOIL'][0, :, :] lat1 = ncid1.variables['lat'][:] lon1 = ncid1.variables['lon'][:] qsoil = -1 * 2.835 * 1e6 * qsoil # use gris as a mask to mask fsa array gris_mask = ncid1.variables['gris_mask'][0, :, :] gris_mask = ma.masked_equal(gris_mask, 0) qsoil_mask = ma.masked_array(qsoil, mask=gris_mask.mask) qsoil = qsoil_mask # read f_racmo_latf_jja, RACMO variables ncid2 = Dataset(f_racmo_latf_jja) latf = ncid2.variables['latf'][0, :, :] # read f_racmo_mask, the lat/lon for RACMO data ncid4 = Dataset(f_racmo_mask) lat4 = ncid4.variables['lat'][:] lon4 = ncid4.variables['lon'][:] racmo_elev = ncid4.variables['Topography'][0, 0, :, :] # Use gris as a mask to mask fsa array. gris_mask = ncid4.variables['GrIS_mask'][:] gris_mask = ma.masked_equal(gris_mask, 0) latf_mask = ma.masked_array(latf, mask=gris_mask.mask) latf = latf_mask # read f_racmo_latf_jja_remapped, the remapped RACMO file to calculate difference ncid3 = Dataset(f_racmo_latf_jja_remapped) remap_latf = ncid3.variables['latf'][0, :, :] diff = qsoil - remap_latf # print(np.max(diff)) # print(np.min(diff)) # ------- PLOT -------- # Open a workstation for drawing the plots wkres = Ngl.Resources() # wkres.wkOrientation = "portrait" # "portrait" or "landscape" wks_type = "png" wks_img = str(os.path.join(out_path, "CESM_RACMO23_latf_JJA")) 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 latf contour of CESM ------- res1 = Ngl.Resources() res1.cnFillPalette = "WhiteBlue" 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.cnLevelSelectionMode = "ExplicitLevels" res1.cnLevels = np.arange(-30, 0, 5) res1.lbLabelBarOn = True # Turn on labelbar. res1.lbLabelFontHeightF = 0.04 # res1.pmLabelBarOrthogonalPosF = -0.05 res1.sfXArray = lon1 res1.sfYArray = lat1 # --- for the latf contour of RACMO ------- res2 = Ngl.Resources() res2.cnFillPalette = "WhiteBlue" 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.cnLevelSelectionMode = "ExplicitLevels" res2.cnLevels = np.arange(-30, 0, 5) res2.lbLabelBarOn = True # Turn on labelbar. res2.lbOrientation = "Vertical" # Verticle labelbar res2.lbLabelFontHeightF = 0.04 # Make fonts smaller. res2.sfXArray = lon4 res2.sfYArray = lat4 # --- for the latf contour of CESM-RACMO ------- res3 = Ngl.Resources() res3.cnFillPalette = "BlueWhiteOrangeRed" 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.cnLevelSelectionMode = "ExplicitLevels" res3.cnLevels = np.arange(-30, 30, 5) res3.lbLabelBarOn = True # Turn on labelbar. res3.lbOrientation = "Vertical" # Verticle labelbar res3.lbLabelFontHeightF = 0.03 # Make fonts smaller. res3.sfXArray = lon1 res3.sfYArray = lat1 # ---- 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 = lon sres.sfYArray = lat # ---- for the elevation of RACMO ------- sres1 = Ngl.Resources() sres1.nglDraw = False # Don't draw individual plots sres1.nglFrame = False # Don't advance frame. sres1.cnFillOn = False sres1.cnLinesOn = True sres1.cnLineLabelsOn = False sres1.cnLevelSelectionMode = "ExplicitLevels" sres1.cnLevels = [0, 1000, 2000, 3000] sres1.sfXArray = lon4 sres1.sfYArray = lat4 # ---- Overlay plots, each one has its own ID # overlay ice on base map, and then overlay elevation on ice # usrf is elevation of 5km CESM, racmo_elev is elevation of RACMO usrf_plot1 = Ngl.contour(wks, usrf, sres) usrf_plot2 = Ngl.contour(wks, racmo_elev, sres1) usrf_plot3 = Ngl.contour(wks, usrf, sres) qsoil_plot = Ngl.contour(wks, qsoil, res1) latf_plot = Ngl.contour(wks, latf, res2) diff_plot = Ngl.contour(wks, diff, res3) # Creat multiple figures and draw, which now contains the elevation and radiation # "[1,3]" indicates 1 row, 3 columns. map_title = ["CESM, LHF", "RACMO, LHF", "CESM-RACMO, LHF"] 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], qsoil_plot) Ngl.overlay(plot[0], usrf_plot1) Ngl.overlay(plot[1], latf_plot) Ngl.overlay(plot[1], usrf_plot2) Ngl.overlay(plot[2], diff_plot) Ngl.overlay(plot[2], usrf_plot3) Ngl.panel(wks, plot, [1, nmap]) 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