# 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 2-meter air temperature (°C) over Greenland for every winter (December-January-Feburary; DJF) 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 = "Winter 2-meter air temperature" 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_model_lnd_climo_djf = os.path.join(config['model_lnd_climos'], 'b.e10.BG20TRCN.f09_g16.002_DJF_climo.nc') f_model_atm_climo_djf = os.path.join(config['model_atm_climos'], 'b.e10.BG20TRCN.f09_g16.002_DJF_climo.nc') f_racmo_t2m_djf = os.path.join(config['racmo_data'], 'racmo23_GRN_monthly.t2m.1980-1999.DJF.nc') f_racmo_t2m_djf_remapped = os.path.join(config['racmo_data'], 'racmo23_GRN_monthly.t2m.1980-1999.remap2cesm.DJF.nc') f_racmo_mask = os.path.join(config['racmo_data'], '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_model_atm_climo_djf and file2, CESM variable ncid1 = Dataset(f_model_atm_climo_djf) w_model = ncid1.variables['TREFHT'][0, :, :] lat1 = ncid1.variables['lat'][:] lon1 = ncid1.variables['lon'][:] # --- original RACMO data without remapping [312,306] ncid2 = Dataset(f_racmo_t2m_djf) 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_djf_remapped) w_racmo = ncid22.variables['t2m'][0, :, :] lat22 = ncid22.variables['lat'][:] lon22 = ncid22.variables['lon'][:] w_model = w_model - 273.15 # convert to celcius w_racmo = w_racmo - 273.15 w_racmo_l = w_racmo_l - 273.15 diff = w_model - w_racmo # locate the greenland area using the GreenLand mask ncid8 = Dataset(f_model_lnd_climo_djf) gris_mask = ncid8.variables['gris_mask'][0, :, :] gris_mask = ma.masked_equal(gris_mask, 0) # mask out non RACMO regions, icemask var in f_model_lnd_climo_djf 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.wkOrientation = "portrait" # "portrait" or "landscape" wks_type = "png" wks_img = str(os.path.join(out_path, "CESM_RACMO23_t2m_DJF")) 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 = "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.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.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.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.cnFillPalette = "BlueWhiteOrangeRed" 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_model, 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(title, ' '.join(describe.split()), img_link) img_elem['Height'] = config['image_height'] img_list.append(img_elem) return img_list