import numpy as np import numpy.ma as ma import Nio, Ngl, os # ---------------- Data source in TITAN ------------------------ input_file0 = '/lustre/atlas1/cli115/world-shared/4ue/Greenland_5km_v1.1_SacksRev_c110629.nc' input_file8 = '/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/lnd/climos/20180215.DECKv1b_H1.ne30_oEC.edison_JJA_climo.nc' input_file1 = '/lustre/atlas1/cli115/world-shared/4ue/20180215.DECKv1b_H1.ne30_oEC.edison/postproc/atm/climos/20180215.DECKv1b_H1.ne30_oEC.edison_JJA_climo.nc' input_file2 = '/lustre/atlas1/cli115/world-shared/4ue/obs_data/ISCCP_ANN_climo.nc' input_file3 = '/lustre/atlas1/cli115/world-shared/4ue/obs_data/CLOUDSAT_ANN_climo.nc' input_file22 = '/lustre/atlas1/cli115/world-shared/4ue/racmo23_GRN_monthly/remapped_racmo/racmo23_GRN_monthly.t2m.1980-1999.remap2cesm.JJA.nc' # -------------------------------------------------------------- # read input_file0, the elevation data ncid0 = Nio.open_file(input_file0) usrf = ncid0.variables['usrf'][0,:,:] lat = ncid0.variables['lat'][0,:,:] lon = ncid0.variables['lon'][0,:,:] # read input_file1 and file2, E3SM variable ncid1 = Nio.open_file(input_file1) cesm = ncid1.variables['CLDHGH'][0,:,:] lat1 = ncid1.variables['lat'][:] lon1 = ncid1.variables['lon'][:] cesm = cesm * 100 # to get percent #--- ISCCP data without remapping [64,128] ncid2 = Nio.open_file(input_file2) isccp_hgh = ncid2.variables['CLDHGH'][0,:,:] lat2 = ncid2.variables['lat'][:] lon2 = ncid2.variables['lon'][:] #--- CLOUDSAT data without remapping [94,192] ncid3 = Nio.open_file(input_file3) cloudsat_hgh = ncid3.variables['CLDHGH'][:,:] lat3 = ncid3.variables['lat'][:] lon3 = ncid3.variables['lon'][:] # locate the greenland area using the GreenLand mask #ncid8 = Nio.open_file(input_file8) #GrisMask = ncid8.variables['gris_mask'][0,:,:] #Gris_mask = ma.masked_equal(GrisMask,0) print "contour plot of CLDHGH vs obs" maxmodel=np.max(cesm) maxisccp=np.max(isccp_hgh) maxcldsat=np.max(cloudsat_hgh) print "Max CLDHGH model",maxmodel print "Max CLDHGH isccp",maxisccp print "Max CLDHGH cldsat",maxcldsat ##------- PLOT -------- # Open a workstation for drawing the plots wkres = Ngl.Resources() #wkres.wkColorMap = "WhiteBlueGreenYellowRed" wkres.wkColorMap = "BlueWhiteOrangeRed" wkres.wkOrientation = "portrait" # "portrait" or "landscape" wks_type = "ps" wks = Ngl.open_wks(wks_type,"/lustre/atlas1/cli115/world-shared/4ue/20171122.beta3rc10_1850.ne30_oECv3_ICG.edison/postproc/atm/plots/E3SM_cldhgh_ANN",wkres) #wks = Ngl.open_wks(wks_type,"E3SM_cldhgh_ANN",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 E3SM 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(0,100,10) res1.lbLabelBarOn = True # Turn on labelbar. res1.lbLabelFontHeightF = 0.04 res1.sfXArray = lon1 res1.sfYArray = lat1 #--- for the data 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(0,100,10) res2.lbLabelBarOn = True # Turn on labelbar. res2.lbOrientation = "Vertical" # Verticle labelbar res2.lbLabelFontHeightF = 0.04 # Make fonts smaller. res2.sfXArray = lon2 res2.sfYArray = lat2 #--- for the data contour ------- 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.trGridType = "TriangularMesh" res3.cnLevelSelectionMode = "ExplicitLevels" res3.cnLevels = np.arange(0,100,10) res3.lbLabelBarOn = True # Turn on labelbar. res3.lbOrientation = "Vertical" # Verticle labelbar res3.lbLabelFontHeightF = 0.04 # Make fonts smaller. res3.sfXArray = lon3 res3.sfYArray = lat3 #---- 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 #---- 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,usrf,sres) usrf_plot3 = Ngl.contour(wks,usrf,sres) cesm_plot = Ngl.contour(wks,cesm,res1) isccp_plot = Ngl.contour(wks,isccp_hgh,res2) cldsat_plot = Ngl.contour(wks,cloudsat_hgh,res3) # Creat multiple figures and draw, which now contains the elevation and temperature # "[1,3]" indicates 1 row, 3 columns. map_title = ["E3SM", "ISCCP","CLOUDSAT"] 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],cesm_plot) Ngl.overlay(plot[0],usrf_plot1) Ngl.overlay(plot[1],isccp_plot) Ngl.overlay(plot[1],usrf_plot2) Ngl.overlay(plot[2],cldsat_plot) Ngl.overlay(plot[2],usrf_plot3) Ngl.panel(wks,plot,[1,3]) Ngl.end()