-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathwind-vector-plot.ncl
134 lines (90 loc) · 3.6 KB
/
wind-vector-plot.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
;;; NCL script for creating wind-vector plots of NARCCAP RCM data
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;; Uncomment the following variables to hard-code them here, or pass values in via command-line. Example:
;; ncl -x ufile=\"uas_RCM3_cgcm3_2066010103.nc\" vfile=\"vas_RCM3_cgcm3_2066010103.nc\" outfile=\"wind_vec\" timestep=0 wind-vector-plot.ncl
;ufile = "uas_RCM3_cgcm3_2066010103.nc"
;vfile = "vas_RCM3_cgcm3_2066010103.nc"
;outfile = "wind_vec"
;timestep = 0 ;; defaults to last step in file
;; open files
uin = addfile(ufile, "r")
vin = addfile(vfile, "r")
;; read in coordinate variables first
lat = uin->lat
lon = uin->lon
time = uin->time
nx = dimsizes(uin->xc)
ny = dimsizes(uin->yc)
;; default to last timestep in file for QC convenience
if (.not. isvar("timestep")) then
timestep = dimsizes(time) - 1
end if
;; default title for convenience
if (.not. isvar("title")) then
date = ut_calendar(time(timestep),0)
timestamp = ""+date(0,0)+"/"+date(0,1)+"/"+date(0,2)+" "+sprintf("%02.0f", date(0,3))+":"+sprintf("%02.0f", date(0,4))
title = systemfunc("basename "+ufile)+"~C~"+systemfunc("basename "+vfile)+"~C~"+timestamp
end if
;; read in data.
udata = uin->uas(timestep,:,:)
vdata = vin->vas(timestep,:,:)
;; hook up coordinate variables
udata@lat2d = lat
udata@lon2d = lon
vdata@lat2d = lat
vdata@lon2d = lon
;; read in map projection info
pname = udata@grid_mapping
proj = uin->$pname$
;; set plotting parameters
res = True
;; map projection stuff.
if (lower_case(proj@grid_mapping_name) .eq. "lambert_conformal_conic") then
res@mpProjection = "LambertConformal"
res@mpLambertMeridianF = proj@longitude_of_central_meridian
res@mpLambertParallel1F = proj@standard_parallel(0)
res@mpLambertParallel2F = proj@standard_parallel(1)
end if
if (lower_case(proj@grid_mapping_name) .eq. "transverse_mercator") then
res@mpProjection = "Mercator"
res@mpCenterLatF = proj@latitude_of_projection_origin
res@mpCenterLonF = proj@longitude_of_central_meridian
end if
if (lower_case(pname) .eq. "polar_stereographic") then
res@mpProjection= "Stereographic"
res@mpCenterLonF= proj@straight_vertical_longitude_from_pole - 360
res@mpCenterLatF = proj@latitude_of_projection_origin
end if
if (lower_case(proj@grid_mapping_name) .eq. "rotated_latitude_longitude") then
res@mpProjection = "Mercator"
res@mpCenterLonF = proj@grid_north_pole_longitude - 180
res@mpCenterLatF = 90 - proj@grid_north_pole_latitude
end if
;; /map projection stuff
;; map boundaries
;; expand slighltly to be sure we plot things on the borders
res@mpLimitMode = "Corners"
res@mpLeftCornerLatF = lat(0,0) -1 ;; SW corner
res@mpLeftCornerLonF = lon(0,0) -1 ;; SW corner
res@mpRightCornerLatF = lat(yc|ny-1,xc|nx-1) + 1 ;; NE corner
res@mpRightCornerLonF = lon(yc|ny-1,xc|nx-1) + 1 ;; NE corner
;; general plotting stuff
res@mpFillOn = False
res@gsnMaximize = True
res@vcGlyphStyle = "CurlyVector"
res@vcRefMagnitudeF = 10.0
res@vcRefLengthF = 0.02
res@vcMinDistanceF = 0.02
res@lbLabelAutoStride = True
res@lbBoxLinesOn = False
res@gsnSpreadColors = True
res@tiMainString = title
;; create the plot
wks = gsn_open_wks("ps", outfile)
gsn_define_colormap(wks,"nrl_sirkes") ;; 21 colors + fg/bg
plot = gsn_csm_vector_map(wks, udata, vdata, res)
exit
;; Copyright 2009-2012 Univ. Corp. for Atmos. Research
;; Author: Seth McGinnis, [email protected]