-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathslp.ncl
executable file
·202 lines (172 loc) · 6.68 KB
/
slp.ncl
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
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
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
; Plots MSLP, Precip, and Winds
begin
; First, let's set some variables. This will dictate how our image files
; are named. File naming convention is similar to Dan Leins' GRaDS
; scripts. If you're using Leins' web app to display images, you should
; not have to do much re-configuring.
domain = "d01_"
domname = "pih_"
prod = "mslp"
directory=str_get_cols(domname,0,2)
; Now, let's load some NetCDF files from our WRF. Change 'dir' to
; the location of your WRF output files.
dir = "/wrf/uems/runs/"+directory+"/wrfprd/"
fils = systemfunc("ls "+dir+"wrfout_"+domain+"*")
a = addfiles(fils+".nc","r")
; On to reading the "cities1.txt" or "cities2.txt" file to plot cities on the map.
index = str_get_cols(domain,2,2)
cities = asciiread("cities"+index+".txt",-1,"string")
delim = ","
lat = tofloat(str_get_field(cities,1,delim))
lon = tofloat(str_get_field(cities,2,delim))
city = str_get_field(cities,3,delim)
; Import attributes from the netCDF file or bad things will happen
; in the next section.
locres = True
locres@MAP_PROJ = a[0]@MAP_PROJ
locres@TRUELAT1 = a[0]@TRUELAT1
locres@TRUELAT2 = a[0]@TRUELAT2
locres@STAND_LON = a[0]@STAND_LON
locres@REF_LAT = a[0]@CEN_LAT
locres@REF_LON = a[0]@CEN_LON
locres@KNOWNI = a[0]@$"WEST-EAST_GRID_DIMENSION"$/2
locres@KNOWNJ = a[0]@$"SOUTH-NORTH_GRID_DIMENSION"$/2
locres@DX = a[0]@DX
locres@DY = a[0]@DY
; Convert WRF lat/lon to NCL-friendly coordinates.
loc = wrf_ll_to_ij(lon,lat,locres)
lo = toint(loc(0,:))
la = toint(loc(1,:))
; Define the file format and the image resolution.
type = "png"
type@wkWidth = 1800
type@wkHeight = 1200
; Set some basic resources
; Title
res = True
res@MainTitle = "Pocatello WRF"
pltres = True
pltres@FramePlot = False
; Map resources
mpres = True
mpres@mpDataBaseVersion = "Ncarg4_1"
mpres@mpOutlineBoundarySets = "AllBoundaries"
mpres@mpGeophysicalLineColor = "Black"
mpres@mpNationalLineColor = "Black"
mpres@mpUSStateLineColor = "Black"
mpres@mpGridLineColor = "Black"
mpres@mpLimbLineColor = "Black"
mpres@mpPerimLineColor = "Black"
mpres@mpCountyLineColor = "Brown"
mpres@mpCountyLineDashPattern = 0
mpres@mpCountyLineThicknessF = 0.5
mpres@mpGeophysicalLineThicknessF = 3.0
mpres@mpGridLineThicknessF = 0.0
mpres@mpLimbLineThicknessF = 2.0
mpres@mpNationalLineThicknessF = 3.0
mpres@mpUSStateLineThicknessF = 3.0
; What times and how many time steps are in the data set?
times = wrf_user_getvar(a,"times",-1) ; get all times in the file
ntimes = dimsizes(times) ; number of times in the file
; Start the time loop
do it = 0,ntimes-1
print("Working on time: " + times(it) )
res@TimeLabel = times(it) ; Set Valid time to use on plots
its = sprintf("%02g",it) ; 2-digit forecast hour for file-naming
; Creating the workstation environment for the plot, define file name
; and select the colormap.
wks = gsn_open_wks(type,domain+domname+prod+its+"_syn")
; gsn_merge_colormaps(wks,"cmocean_gray","precip2_17lev")
; gsn_define_colormap(wks,"precip2_17lev")
; What data are we putting on the map?
slp = wrf_user_getvar(a,"slp",it) ; MSLP
hgt = wrf_user_getvar(a,"ter",it) ; Terrain
terrain=hgt*3.281 ; Convert to ft.
rainnc = wrf_user_getvar(a,"RAINNC",it); Precip
rainc = wrf_user_getvar(a,"RAINC",it) ; Some other precip variable
precip = rainnc + rainc
if(it .gt. 0) then
rainnc_old=wrf_user_getvar(a,"RAINNC",it-1)
rainc_old=wrf_user_getvar(a,"RAINC",it-1)
precip := precip-rainnc_old-rainc_old
end if
pcp = precip*0.0393701 ; Convert to inches
u10 = wrf_user_getvar(a,"U10",-1) ; 10m u wind
v10 = wrf_user_getvar(a,"V10",-1) ; 10m v wind
u10 = u10* 1.94384
v10 = v10* 1.94384
; Set resources for these variables.
slp@description = "Sea Level Pressure"
slp@units = "hPa"
terrain@description = "Terrain Height"
terrain@units = "ft"
pcp@description = "Precipitation"
pcp@units = "inches"
u10@description = "10m Wind"
u10@units = "Knots"
cmap = read_colormap_file("precip2_17lev")
cmap(0,:) = (/0,0,0,0/)
; Moving label values closer to the color bar.
res@lbLabelOffsetF = 0.05
; Plotting options for Precip
opts = res
opts@ContourParameters = (/0., 0.5, 0.01/)
opts@SubFieldTitle = " Max: "+decimalPlaces(max(pcp),2,True)
opts@Footer = False
opts@cnFillOn = True
opts@cnFillPalette = cmap
contour_precip=wrf_contour(a[it],wks,pcp,opts)
delete(opts)
; Plotting options for Terrain
opts = res
opts@cnFillOn = True
opts@ContourParameters = (/0., 3000., 100./)
opts@cnFillPalette = "MPL_terrain"
opts@lbLabelsOn = False
opts@lbLabelBarOn = False
opts@Footer = False
contour_terrain=wrf_contour(a[it],wks,terrain,opts)
delete(opts)
; Plotting options for SLP
opts = res
opts@cnLineColor = "Black"
opts@cnLineThicknessF = 3.0
opts@SubFieldTitle = " Max: "+toint(max(slp))+" Min: "+toint(min(slp))
opts@Footer = False
opts@lbBottomMarginF = -0.2
opts@ContourParameters = (/ 1 /)
; opts@gsnSpreadColorEnd = -3 ; End third from the last color in color map
wrf_smooth_2d(slp,40)
contour_slp = wrf_contour(a[it],wks,slp,opts)
delete(opts) ; delete opts before adding another contour.
; Plotting options for Wind
opts = res
opts@Footer = False
opts@vcLevelColors = "yellow"
opts@vcRefLengthF = 0.02
opts@vcWindBarbLineThicknessF = 2.0
vector = wrf_vector(a[it],wks,u10(it,:,:),v10(it,:,:),opts)
delete(opts)
; Time to overlay your plots onto the map.
plot = wrf_map_overlays(a[it],wks,(/contour_terrain, contour_precip, contour_slp, vector/),pltres,mpres)
; Now we can plot cities on the map.
; places = dimsizes(cities) ; How many cities in our cities.txt file?
; text = new(places,float) ; Now we sample some data and put it under
; Now set resources for the text
; txres = True
; txres@txFontHeightF = 0.008
; define an offset based on the domain
; if (domain .eq. "d01_") then
; offset = 0.1
; else
; offset = 0.05
; end if
; gsn_text(wks,plot,city,lon,lat+offset,txres) ; Plot the city with a slight offset.
; gsn_text(wks,plot,".",lon,lat,txres) ; Plots a dot at the city location.
; gsn_text(wks,plot,text,lon,lat-offset,txres) ; Plots the value below
;draw(plot)
frame(wks)
system("convert -trim "+domain+domname+prod+its+"_syn.png "+domain+domname+prod+its+"_syn.png")
system("convert -border 10 -bordercolor white "+domain+domname+prod+its+"_syn.png "+domain+domname+prod+its+"_syn.png")
end do ; END OF TIME LOOP
end