# 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 summer (June-July-August; JJA) average net radiation (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). """ 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_swsn_jja = os.path.join(config['racmo_data'], 'racmo23_GRN_monthly.swsn.1980-1999.JJA.nc') f_racmo_swsn_jja_remapped = os.path.join(config['racmo_data'], 'racmo23_GRN_monthly.swsn.1980-1999.remap2cesm.JJA.nc') f_racmo_lwsn_jja = os.path.join(config['racmo_data'], 'racmo23_GRN_monthly.lwsn.1980-1999.JJA.nc') f_racmo_lwsn_jja_remapped = os.path.join(config['racmo_data'], 'racmo23_GRN_monthly.lwsn.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 [:] 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 varialbes ncid1 = Dataset(f_cesm_lnd_climo_jja) fsa = ncid1.variables['FSA'][0, :, :] lat1 = ncid1.variables['lat'][:] lon1 = ncid1.variables['lon'][:] fira = ncid1.variables['FIRA'][0, :, :] rnet_cesm = fsa - fira # use gris as a mask to mask fsa array gris_mask = ncid1.variables['gris_mask'][0, :, :] gris_mask = ma.masked_equal(gris_mask, 0) rnet_cesm_mask = ma.masked_array(rnet_cesm, mask=gris_mask.mask) rnet_cesm = rnet_cesm_mask # read input_file2, RACMO variables ncid21 = Dataset(f_racmo_swsn_jja) swsn = ncid21.variables['swsn'][0, :, :] ncid22 = Dataset(f_racmo_lwsn_jja) lwsn = ncid22.variables['lwsn'][0, :, :] rnet = swsn + lwsn # 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) rnet_mask = ma.masked_array(rnet, mask=gris_mask.mask) rnet = rnet_mask # read input_file3, the remapped RACMO file to calculate difference ncid31 = Dataset(f_racmo_swsn_jja_remapped) remap_swsn = ncid31.variables['swsn'][0, :, :] ncid32 = Dataset(f_racmo_lwsn_jja_remapped) remap_lwsn = ncid32.variables['lwsn'][0, :, :] remap_rnet = remap_swsn + remap_lwsn diff = rnet_cesm - remap_rnet # 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_rnet_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 rnet 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(-20, 100, 20) res1.lbLabelBarOn = True # Turn on labelbar. res1.lbLabelFontHeightF = 0.04 # res1.pmLabelBarOrthogonalPosF = -0.05 res1.sfXArray = lon1 res1.sfYArray = lat1 # --- for the rnet 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(-20, 100, 20) 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 rnet 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(-80, 40, 20) 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) rnet_cesm_plot = Ngl.contour(wks, rnet_cesm, res1) rnet_plot = Ngl.contour(wks, rnet, 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, R~B3~net", "RACMO, R~B3~net", "CESM-RACMO, R~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], rnet_cesm_plot) Ngl.overlay(plot[0], usrf_plot1) Ngl.overlay(plot[1], rnet_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_rnet', ' '.join(describe.split()), img_link) img_elem['Height'] = config['image_height'] img_list.append(img_elem) return img_list