cesm_cldlow_ann.py 6.92 KB
Newer Older
1
2
3
4
5
6
import os
import Ngl
import numpy as np

from netCDF4 import Dataset

7
from livvkit.util import elements as el
8

9
10
11
12
13

describe = """CESM_cldlow_ANN plot."""


def make_plot(config=None, out_path='.',
14
15
16
17
18
19
20
21
22
23
              cloud_path='/lustre/atlas1/cli115/world-shared/4ue/obs_data/',
              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_atm_climo_jja = os.path.join(cesm_path, 'postproc/atm/climos/b.e10.BG20TRCN.f09_g16.002_JJA_climo.nc')
    f_isccp = os.path.join(cloud_path, 'ISCCP_ANN_climo.nc')
    f_cloudsat = os.path.join(cloud_path, 'CLOUDSAT_ANN_climo.nc')
    # --------------------------------------------------------------

24
25
    img_list = []

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
    # read f_cism, the elevation data
    cism = Dataset(f_cism)
    cism_usrf = cism.variables['usrf'][0, :, :]
    cism_lat = cism.variables['lat'][0, :, :]
    cism_lon = cism.variables['lon'][0, :, :]

    # read f_cesm_atm_climo_jja, CESM variable
    cesm = Dataset(f_cesm_atm_climo_jja)
    cesm_low = cesm.variables['CLDLOW'][0, :, :]
    cesm_lat = cesm.variables['lat'][:]
    cesm_lon = cesm.variables['lon'][:]

    cesm_low = cesm_low * 100  # to get percent

    # --- ISCCP data without remapping [64,128]
    isccp = Dataset(f_isccp)
    isccp_low = isccp.variables['CLDLOW'][0, :, :]
    isccp_lat = isccp.variables['lat'][:]
    isccp_lon = isccp.variables['lon'][:]

    # --- CLOUDSAT data without remapping [94,192]
    cloudsat = Dataset(f_cloudsat)
    cloudsat_low = cloudsat.variables['CLDLOW'][:, :]
    cloudsat_lat = cloudsat.variables['lat'][:]
    cloudsat_lon = cloudsat.variables['lon'][:]

52
53
54
55
56
57
58
    # print("contour plot of CLDLOW vs obs")
    # maxmodel = np.max(cesm_low)
    # maxisccp = np.max(isccp_low)
    # maxcldsat = np.max(cloudsat_low)
    # print("Max CLDLOW model:  {}".format(maxmodel))
    # print("Max CLDLOW isccp:  {}".format(maxisccp))
    # print("Max CLDLOW cldsat: {}".format(maxcldsat))
59
60
61
62
63
64

    # ------- PLOT --------
    #  Open a workstation for drawing the plots
    wkres = Ngl.Resources()
    # wkres.wkColorMap = "WhiteBlueGreenYellowRed"
    wkres.wkColorMap = "BlueWhiteOrangeRed"
65
    # wkres.wkOrientation = "portrait"  # "portrait" or "landscape"
66
    wks_type = "png"
67
    wks_img = str(os.path.join(out_path, "CESM_cldlow_ANN"))
68
    wks = Ngl.open_wks(wks_type, wks_img, wkres)
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
185
186
187

    # --- 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(0, 100, 10)
    res1.lbLabelBarOn = True  # Turn on labelbar.
    res1.lbLabelFontHeightF = 0.04

    res1.sfXArray = cesm_lon
    res1.sfYArray = cesm_lat

    # --- 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 = isccp_lon
    res2.sfYArray = isccp_lat

    # --- 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 = cloudsat_lon
    res3.sfYArray = cloudsat_lat

    # ---- 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 = cism_lon
    sres.sfYArray = cism_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, cism_usrf, sres)
    usrf_plot2 = Ngl.contour(wks, cism_usrf, sres)
    usrf_plot3 = Ngl.contour(wks, cism_usrf, sres)

    cesm_plot = Ngl.contour(wks, cesm_low, res1)
    isccp_plot = Ngl.contour(wks, isccp_low, res2)
    cldsat_plot = Ngl.contour(wks, cloudsat_low, res3)

    # Creat multiple figures and draw, which now contains the elevation and temperature
    # "[1,3]" indicates 1 row, 3 columns.
    map_title = ["CESM", "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])

188
189
190
191
192
193
194
195
196
197
198
    img_link = os.path.join(os.path.basename(out_path),
                            os.path.basename(wks_img + '.' + wks_type))
    img_elem = el.image('CESM_cldlow_ANN',
                        ' '.join(describe.split()),
                        img_link)
    if config:
        img_elem['Height'] = config['image_height']
    img_list.append(img_elem)

    return img_list

199
200
201

if __name__ == '__main__':
    make_plot()