model_cldtot_ann.py 6.78 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
describe = """
Climatological annual average of cloud cover over Greenland 
12
for Model (left), ISCCP (middle; Rossow and Schiffer, 1999), and CLOUDSAT (right; 
13
Kay and Gettelman, 2009).  The black solid lines denote the 0, 1000, 2000, and 
14
3000 meter Greenland ice sheet elevation contours as seen by the Model.
15
"""
16

17
18
title = "Cloud cover over Greenland"

19

20
21
22
def make_plot(config, out_path='.'):
    f_isccp = os.path.join(config['cloud_data'], 'ISCCP_ANN_climo.nc')
    f_cloudsat = os.path.join(config['cloud_data'], 'CLOUDSAT_ANN_climo.nc')
23
24
25

    img_list = []

26
    cism = Dataset(config['glc_surf'])
27
28
29
30
31
    cism_usrf = cism.variables['usrf'][0, :, :]
    cism_lat = cism.variables['lat'][0, :, :]
    cism_lon = cism.variables['lon'][0, :, :]
    
    # read f_cesm_atm_climo_jja and file2, CESM variable
32
33
34
35
    model = Dataset(config['atm_climo'])
    model_tot = model.variables['CLDTOT'][0, :, :]
    model_lat = model.variables['lat'][:]
    model_lon = model.variables['lon'][:]
36
    
37
    model_tot = model_tot * 100  # to get percent
38
39
40
41
42
43
44
45
46
47
48
49
50
    
    # --- ISCCP data without remapping [64,128]
    isccp = Dataset(f_isccp)
    isccp_tot = isccp.variables['CLDTOT'][0, :, :]
    isccp_lat = isccp.variables['lat'][:]
    isccp_lon = isccp.variables['lon'][:]
    
    # --- CLOUDSAT data without remapping [94,192]
    cloudsat = Dataset(f_cloudsat)
    cloudsat_tot = cloudsat.variables['CLDTOT'][:, :]
    cloudsat_lat = cloudsat.variables['lat'][:]
    cloudsat_lon = cloudsat.variables['lon'][:]
    
51
    # print("contour plot of CLDTOT vs obs")
52
    # maxmodel = np.max(model_tot)
53
54
55
56
57
    # maxisccp = np.max(isccp_tot)
    # maxcldsat = np.max(cloudsat_tot)
    # print("Max CLDTOT model:  {}".format(maxmodel))
    # print("Max CLDTOT isccp:  {}".format(maxisccp))
    # print("Max CLDTOT cldsat: {}".format(maxcldsat))
58
59
60
    
    # ------- PLOT --------
    #  Open a workstation for drawing the plots
61
    cmap = ["grey98","grey94","grey90","grey86","grey80","grey74","grey68","grey60","grey52","grey44","grey36"]
62
    wkres = Ngl.Resources()
63
    # wkres.wkOrientation = "portrait"  # "portrait" or "landscape"
64
    wks_type = "png"
65
    wks_img = str(os.path.join(out_path, "Model_cldtot_ANN"))
66
    wks = Ngl.open_wks(wks_type, wks_img, wkres)
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
    
    # --- 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
    
88
    # --- for the Model contour -------
89
    res1 = Ngl.Resources()
90
    res1.cnFillPalette = "percent_11lev"
91
92
93
94
95
96
97
98
99
100
101
102
    res1.nglDraw = False  # Don't draw individual plots
    res1.nglFrame = False  # Don't advance frame.
    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
    
103
104
    res1.sfXArray = model_lon
    res1.sfYArray = model_lat
105
106
107
    
    # --- for the data contour -------
    res2 = Ngl.Resources()
108
    res2.cnFillPalette = "percent_11lev"
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
    res2.nglDraw = False  # Don't draw individual plots
    res2.nglFrame = False  # Don't advance frame.
    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()
127
    res3.cnFillPalette = "percent_11lev"
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
    res3.nglDraw = False  # Don't draw individual plots
    res3.nglFrame = False  # Don't advance frame.
    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
151
    sres.trGridType = "TriangularMesh"
152
153
154
155
156
157
158
159
160
161
162
    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)
    
163
    model_plot = Ngl.contour(wks, model_tot, res1)
164
165
166
167
168
    isccp_plot = Ngl.contour(wks, isccp_tot, res2)
    cldsat_plot = Ngl.contour(wks, cloudsat_tot, res3)
    
    # Creat multiple figures and draw, which now contains the elevation and temperature
    # "[1,3]" indicates 1 row, 3 columns.
169
    map_title = ["Model", "ISCCP", "CLOUDSAT"]
170
171
172
173
174
175
176
177
    
    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.
178
    Ngl.overlay(plot[0], model_plot)
179
180
181
182
183
184
185
186
    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])

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

    return img_list