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_t2m_ann 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_ann = os.path.join(cesm_path, 'postproc/lnd/climos/b.e10.BG20TRCN.f09_g16.002_ANN_climo.nc') f_cesm_atm_climo_ann = os.path.join(cesm_path, 'postproc/atm/climos/b.e10.BG20TRCN.f09_g16.002_ANN_climo.nc') f_racmo_t2m_ann = os.path.join(racmo_path, 'climos/racmo23_GRN_monthly.t2m.1980-1999.ANN.nc') f_racmo_t2m_ann_remapped = os.path.join(racmo_path, 'remapped_racmo/racmo23_GRN_monthly.t2m.1980-1999.remap2cesm.ANN.nc') f_racmo_mask = os.path.join(racmo_path, 'RACMO23_masks_ZGRN11.nc') # -------------------------------------------------------------- img_list = [] # read f_cism, the elevation data ncid0 = Dataset(f_cism) usrf = ncid0.variables['usrf'][0, :, :] lat = ncid0.variables['lat'][0, :, :] lon = ncid0.variables['lon'][0, :, :] # read f_cesm_atm_climo_ann and file2, CESM variable ncid1 = Dataset(f_cesm_atm_climo_ann) w_cesm = ncid1.variables['TREFHT'][0, :, :] lat1 = ncid1.variables['lat'][:] lon1 = ncid1.variables['lon'][:] # --- original RACMO data without remapping [312,306] ncid2 = Dataset(f_racmo_t2m_ann) w_racmo_l = ncid2.variables['t2m'][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, :, :] # RACMO data after remapping [192,288] ncid22 = Dataset(f_racmo_t2m_ann_remapped) w_racmo = ncid22.variables['t2m'][0, :, :] lat22 = ncid22.variables['lat'][:] lon22 = ncid22.variables['lon'][:] w_cesm = w_cesm - 273.15 # convert to celcius w_racmo = w_racmo - 273.15 w_racmo_l = w_racmo_l - 273.15 diff = w_cesm - w_racmo # locate the greenland area using the GreenLand mask ncid8 = Dataset(f_cesm_lnd_climo_ann) gris_mask = ncid8.variables['gris_mask'][0, :, :] gris_mask = ma.masked_equal(gris_mask, 0) # mask out non RACMO regions, icemask var in f_cesm_lnd_climo_ann diff_mask = ma.masked_array(diff, mask=gris_mask.mask) # print np.max(diff_mask) # print np.min(diff_mask) # ------- PLOT -------- # Open a workstation for drawing the plots wkres = Ngl.Resources() wkres.wkColorMap = "BlueWhiteOrangeRed" # wkres.wkOrientation = "portrait" # "portrait" or "landscape" wks_type = "png" wks_img = str(os.path.join(out_path, "CESM_RACMO23_t2m_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(-44, 12, 4) res1.lbLabelBarOn = True # Turn on labelbar. res1.lbLabelFontHeightF = 0.04 res1.sfXArray = lon1 res1.sfYArray = lat1 # --- for the RACMO 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(-44, 12, 4) 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 diff=CESM-Remapped_RACMO contour ------- res22 = Ngl.Resources() res22.nglDraw = False # Don't draw individual plots res22.nglFrame = False # Don't advance frame. res22.cnLineLabelsOn = False res22.cnFillOn = True res22.cnLinesOn = False res22.cnLineLabelsOn = False res22.cnFillMode = "RasterFill" res22.trGridType = "TriangularMesh" res22.cnLevelSelectionMode = "ExplicitLevels" res22.cnLevels = np.arange(-12, 12, 2) res22.lbLabelBarOn = True # Turn on labelbar. res22.lbOrientation = "Vertical" # Verticle labelbar res22.lbLabelFontHeightF = 0.04 # Make fonts smaller. res22.sfXArray = lon22 res22.sfYArray = lat22 # ---- 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_plot1 = Ngl.contour(wks, usrf, sres) usrf_plot2 = Ngl.contour(wks, racmo_elev, sres1) usrf_plot3 = Ngl.contour(wks, usrf, sres) t2mw_plot = Ngl.contour(wks, w_cesm, res1) t2ms_l_plot = Ngl.contour(wks, w_racmo_l, res2) diff_plot = Ngl.contour(wks, diff_mask, res22) # diff_plot = Ngl.contour(wks,w_racmo,res2) #Remapped RACMO # Creat multiple figures and draw, which now contains the elevation and temperature # "[1,3]" indicates 1 row, 3 columns. map_title = ["CESM T(~S1~o C )", "RACMO T(~S1~o C )", "CESM-RACMO T(~S1~o C )"] 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], t2mw_plot) Ngl.overlay(plot[0], usrf_plot1) Ngl.overlay(plot[1], t2ms_l_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_t2m_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()