E3SM_cldhgh_ANN.py 6.52 KB
Newer Older
1
2
3
4
5
6
7
8

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'
9
10
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'
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
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()