-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathTRANS_slice.py
129 lines (108 loc) · 3.89 KB
/
TRANS_slice.py
1
2
3
4
5
6
7
8
9
10
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
#
# File:
# TRANS_slice.py
#
# Synopsis:
# Illustrates how to create a slice plot
#
# Categories:
# contour plot
#
# Author:
# Karin Meier-Fleischer, based on NCL example
#
# Date of initial publication:
# September 2018
#
# Description:
# This example shows how to create a slice plot.
# xarray is used to read the NetCDF file, but PyNIO can also
# be used. See the lines commented with ###.
#
# Effects illustrated:
# o Read netCDF data
# o Drawing a slice plot
#
# Output:
# A single visualization is produced.
#
# Notes: The data for this example can be downloaded from
# http://www.ncl.ucar.edu/Document/Manuals/NCL_User_Guide/Data/
#
'''
Transition Guide PyNGL Example: TRANS_slice.py
- Read netCDF data
- Drawing a slice plot
18-09-04 kmf
'''
import numpy as np
import Ngl
import xarray as xr
###import Nio
import os
#-- define variables
fname = "rectilinear_grid_3D.nc" #-- data file name
#-- open file and read variables
f = xr.open_dataset(fname) #-- open data file
t = f.t[0,:,::-1,:] #-- first time step, reverse latitude
lev = f.lev.values #-- all levels
lat = f.lat[::-1] #-- reverse latitudes
lon = f.lon.values #-- all longitudes
###f = Nio.open_file(fname,"r") #-- open data file
###t = f.variables["t"][0,:,::-1,:] #-- first time step, reverse latitude
###lev = f.variables["lev"][:] #-- all levels
###lat = f.variables["lat"][::-1] #-- reverse latitudes
###lon = f.variables["lon"][:] #-- all longitudes
nlat = len(lat) #-- number of latitudes
longname = t.long_name
units = f'deg{t.units}'
###longname = f.variables["t"].attributes['long_name']
###units = f.variables["t"].attributes['units']
# Subselect the data around latitude 40 and 41.
t40 = t.sel(lat=slice(40,41)).squeeze()
print(t40.shape) # 17 x 192
#-- open a workstation
wks = Ngl.open_wks("png",os.path.basename(__file__).split('.')[0])
#-- set resources
res = Ngl.Resources() #-- generate an res object for plot
res.nglFrame = False
#-- viewport resources
res.vpXF = 0.1 #-- start x-position of viewport
res.vpYF = 0.9 #-- start y-position of viewport
res.vpWidthF = 0.7 #-- width of viewport
res.vpHeightF = 0.6 #-- height of viewport
#-- contour resources
res.cnFillOn = True #-- turn on contour fill
res.cnFillPalette = "temp_diff_18lev"
res.cnLineLabelsOn = False #-- turn off line labels
res.cnInfoLabelOn = False #-- turn off info label
res.cnLevelSelectionMode = "ManualLevels"#-- select manual levels
res.cnMinLevelValF = 200. #-- minimum contour value
res.cnMaxLevelValF = 290. #-- maximum contour value
res.cnLevelSpacingF = 5. #-- contour increment
res.tiMainString = 'slice at lat = {:4.2f}'.format(t40.lat.values)
res.tiYAxisString = f"{longname} [hPa]"
#-- grid resources
res.sfXArray = lon #-- scalar field x
res.sfYArray = lev/100 #-- scalar field y
#-- reverse y-axis
res.trYReverse = True #-- reverse the Y axis
res.nglYAxisType = "LogAxis" #-- y axis log
#-- draw slice contour plot
plot = Ngl.contour(wks,t40,res)
#-- Retrieve some resources from map for adding labels
vpx = Ngl.get_float(plot,'vpXF')
vpy = Ngl.get_float(plot,'vpYF')
vpw = Ngl.get_float(plot,'vpWidthF')
fnth = Ngl.get_float(plot,'tmXBLabelFontHeightF')
#-- write variable long_name and units to the plot
txres = Ngl.Resources()
txres.txFontHeightF = fnth
txres.txJust = "CenterLeft"
Ngl.text_ndc(wks,longname,vpx,vpy+0.02,txres)
txres.txJust = "CenterRight"
Ngl.text_ndc(wks,units, vpx+vpw,vpy+0.02,txres)
#-- advance the frame
Ngl.frame(wks)
#-- done
Ngl.end()