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 = """CESM_RACMO23_lwsn plot.""" def make_plot(config=None, out_path='.', racmo_path='/lustre/atlas1/cli115/world-shared/4ue/racmo23_GRN_monthly/', 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_lnd_climo_jja = os.path.join(cesm_path, 'postproc/lnd/climos/b.e10.BG20TRCN.f09_g16.002_JJA_196006_200508_climo.nc') f_racmo_lwsn_jja = os.path.join(racmo_path, 'climos/racmo23_GRN_monthly.lwsn.1980-1999.JJA.nc') f_racmo_lwsn_jja_remapped = os.path.join(racmo_path, 'remapped_racmo/racmo23_GRN_monthly.lwsn.1980-1999.remap2cesm.JJA.nc') f_racmo_mask = os.path.join(racmo_path, 'RACMO23_masks_ZGRN11.nc') img_list = [] # -------------------------------------------------------------- # f_cism get the following, matrixsize[301,561] # usrf(0,:,:), lat(0,:,:), lon(0,:,:) # # f_cesm_lnd_climo_jja get FSDS (mean downwelling solar flux),FSA (mean net solar flux, absorbed) # FSR (the mean net solar flux, reflected) # FSDS(0,:,:), FSA(0,:,:), FSR(0,:,:) # Goal: plot 3 figures, radiation overlaying on elevation usrf # read f_cism, elevation of 5km CESM # adding the [:] transform 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) fira = ncid1.variables['FIRA'][0, :, :] lat1 = ncid1.variables['lat'][:] lon1 = ncid1.variables['lon'][:] fira = -1 * fira # use gris as a mask to mask fira array gris_mask = ncid1.variables['gris_mask'][0, :, :] gris_mask = ma.masked_equal(gris_mask, 0) fira_mask = ma.masked_array(fira, mask=gris_mask.mask) fira = fira_mask # read f_racmo_lwsn_jja, RACMO variables ncid2 = Dataset(f_racmo_lwsn_jja) lwsn = ncid2.variables['lwsn'][0, :, :] # lat2 = ncid2.variables['lat'][:] # lon2 = ncid2.variables['lon'][:] # 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 fira array. gris_mask = ncid4.variables['GrIS_mask'][:] gris_mask = ma.masked_equal(gris_mask, 0) lwsn_mask = ma.masked_array(lwsn, mask=gris_mask.mask) lwsn = lwsn_mask # read f_racmo_lwsn_jja_remapped, the remapped RACMO file to calculate difference ncid3 = Dataset(f_racmo_lwsn_jja_remapped) remap_lwsn = ncid3.variables['lwsn'][0, :, :] diff = fira - remap_lwsn # print(np.max(diff)) # print(np.min(diff)) # ------- PLOT -------- # Open a workstation for drawing the plots wkres = Ngl.Resources() # wkres.wkColorMap = "matlab_jet" wkres.wkColorMap = "BlueWhiteOrangeRed" # wkres.wkOrientation = "portrait" # "portrait" or "landscape" wks_type = "png" wks_img = str(os.path.join(out_path, "CESM_RACMO23_lwsn_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 lwsn contour of CESM ------- 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.cnLevelSelectionMode = "ExplicitLevels" res1.cnLevels = np.arange(-70, -20, 5) res1.lbLabelBarOn = True # Turn on labelbar. res1.lbLabelFontHeightF = 0.04 # res1.pmLabelBarOrthogonalPosF = -0.05 res1.sfXArray = lon1 res1.sfYArray = lat1 # --- for the lwsn contour of RACMO ------- 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.cnLevelSelectionMode = "ExplicitLevels" res2.cnLevels = np.arange(-70, -20, 5) res2.lbLabelBarOn = True # Turn on labelbar. res2.lbOrientation = "Vertical" # Verticle labelbar # res2.pmLabelBarHeightF = 0.5 # Change height of labelbar # res2.pmLabelBarWidthF = 0.2 # Change width of labelbar # res2.pmLabelBarOrthogonalPosF = -0.02 # Move labelbar closer to plot res2.lbLabelFontHeightF = 0.04 # Make fonts smaller. res2.sfXArray = lon4 res2.sfYArray = lat4 # --- for the lwsn contour of CESM-RACMO ------- 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.cnLevelSelectionMode = "ExplicitLevels" res3.cnLevels = np.arange(-30, 40, 10) res3.lbLabelBarOn = True # Turn on labelbar. res3.lbOrientation = "Vertical" # Verticle labelbar res3.lbLabelFontHeightF = 0.04 # 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) fira_plot = Ngl.contour(wks, fira, res1) lwsn_plot = Ngl.contour(wks, lwsn, 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, LW~B3~net", "RACMO, LW~B3~net", "CESM-RACMO, LW~B3~net"] 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], fira_plot) Ngl.overlay(plot[0], usrf_plot1) Ngl.overlay(plot[1], lwsn_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('CESM_RACMO23_lwsn', ' '.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()