| 3 |
michaesp |
1 |
;-------------------------------------------------------------------------------------;
|
|
|
2 |
; Trajectory Plot
|
|
|
3 |
; Assumes Lagranto lsl-file style output with five
|
|
|
4 |
; header lines
|
|
|
5 |
;-------------------------------------------------------------------------------------;
|
|
|
6 |
|
|
|
7 |
;-------------------------------------------------------------------------------------;
|
|
|
8 |
; Loading Libraries
|
|
|
9 |
;-------------------------------------------------------------------------------------;
|
|
|
10 |
|
|
|
11 |
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
|
|
|
12 |
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
|
|
|
13 |
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
|
|
|
14 |
|
|
|
15 |
;-------------------------------------------------------------------------------------;
|
|
|
16 |
; Read trajectory file
|
|
|
17 |
;-------------------------------------------------------------------------------------;
|
|
|
18 |
|
|
|
19 |
; --------- Name of trajectory file --------------------------------------------------;
|
|
|
20 |
filename="trajectory.lsl"
|
|
|
21 |
|
|
|
22 |
;---------- Read file ----------------------------------------------------------------;
|
|
|
23 |
delim = " " ; Delimiting character
|
|
|
24 |
|
|
|
25 |
nhead = 5 ; # header lines
|
|
|
26 |
head = readAsciiHead(filename,nhead-1) ; Read header lines
|
|
|
27 |
ncol = str_fields_count(head(2),delim) ; # columns
|
|
|
28 |
|
|
|
29 |
data = asciiread(filename,-1,"float") ; Read data in as float
|
|
|
30 |
n = dimsizes(data); ; Number of entries in data
|
|
|
31 |
ifirst = 4 ; First data element
|
|
|
32 |
ifinal = n-1 ; Last data element
|
|
|
33 |
|
|
|
34 |
timestp = data(ifirst+ncol) - data(ifirst) ; Time step of trajectories
|
|
|
35 |
period = data(ifinal-ncol+1) ; Time period of trajectories
|
|
|
36 |
ntim = floattoint( period/timestp + 1 ) ; # times
|
|
|
37 |
ntra = floattoint( (n-ifirst) / (ntim * ncol ) ) ; # trajectories
|
|
|
38 |
|
|
|
39 |
traj = onedtond(data(ifirst:ifinal),(/ntra,ntim,ncol/)) ; 3d trajectory array
|
|
|
40 |
vars = str_split(head(2),delim) ; Field names (columns)
|
|
|
41 |
|
|
|
42 |
;-------------- Add meta information --------------------------------------------------;
|
|
|
43 |
|
|
|
44 |
; Add some attributes to the variable
|
|
|
45 |
traj@HOUR = stringtoint( str_get_cols(head(0),22,23) )
|
|
|
46 |
traj@DAY = stringtoint( str_get_cols(head(0),19,20) )
|
|
|
47 |
traj@MONTH = stringtoint( str_get_cols(head(0),17,18) )
|
|
|
48 |
traj@YEAR = stringtoint( str_get_cols(head(0),13,16) )
|
|
|
49 |
traj@FIELD_NAMES = vars
|
|
|
50 |
traj@MODEL = "Lagranto"
|
|
|
51 |
|
|
|
52 |
; Set names of dimensions and assign values
|
|
|
53 |
traj!0 = "TRA"
|
|
|
54 |
traj!1 = "TIME"
|
|
|
55 |
traj!2 = "FIELD"
|
|
|
56 |
traj&TIME = fspan(0,period,ntim)
|
|
|
57 |
|
|
|
58 |
;-------------------------------------------------------------------------------------;
|
|
|
59 |
; Plot trajectories
|
|
|
60 |
;-------------------------------------------------------------------------------------;
|
|
|
61 |
|
|
|
62 |
; --------------- Open Workstation and Colormap --------------------------------------
|
|
|
63 |
|
|
|
64 |
wks = gsn_open_wks("png","tra") ; open workstation
|
|
|
65 |
gsn_define_colormap(wks,"temp_19lev") ; use the BlueDarkRed18 colormap
|
|
|
66 |
cmap = gsn_retrieve_colormap(wks) ; colormap
|
|
|
67 |
|
|
|
68 |
; --------------- General resources for plot -----------------------------------------
|
|
|
69 |
|
|
|
70 |
mres = True ; map resource
|
|
|
71 |
mres@gsnDraw = False ; don't draw
|
|
|
72 |
mres@gsnFrame = False ; don't advance frame
|
|
|
73 |
mres@tiMainString = "Trajectories" ; set the main title
|
|
|
74 |
mres@tiMainFont = "0_times_roman" ; font of main title
|
|
|
75 |
mres@gsnMaximize = True ; Maximize plot
|
|
|
76 |
|
|
|
77 |
|
|
|
78 |
; ----------------- World map --------------------------------------------------------
|
|
|
79 |
|
|
|
80 |
mres@mpMaxLatF = 70 ; choose subregion
|
|
|
81 |
mres@mpMinLatF = 0 ; of world map
|
|
|
82 |
mres@mpMaxLonF = 100
|
|
|
83 |
mres@mpMinLonF = 0
|
|
|
84 |
mres@mpFillOn = False ; color land
|
|
|
85 |
mres@mpOutlineOn = True ; outline land/sea boarders
|
|
|
86 |
|
|
|
87 |
map = gsn_csm_map_ce(wks,mres) ; Draw world map
|
|
|
88 |
draw(map)
|
|
|
89 |
|
|
|
90 |
; ----------------- Polyline (including coloring of trajectories ) --------------------
|
|
|
91 |
|
|
|
92 |
|
|
|
93 |
; --- Specify the coloring field
|
|
|
94 |
cnField = "p" ; Coloring according to pressure
|
|
|
95 |
cnIndex = 3 ; Column index for coloring field
|
|
|
96 |
cnLevels = fspan(200,1000,17) ; Coloring levels
|
|
|
97 |
cnLabels = flt2string(cnLevels) ; Labels
|
|
|
98 |
|
|
|
99 |
|
|
|
100 |
; --- Resources for the polylines (except for color)
|
|
|
101 |
pres = True ; polyline resource
|
|
|
102 |
pres@tfDoNDCOverlay = True
|
|
|
103 |
pres@gsLineThicknessF = 3.0 ; line thickness
|
|
|
104 |
|
|
|
105 |
|
|
|
106 |
; --- Resources for the colorbar/labelbar
|
|
|
107 |
lres = True
|
|
|
108 |
lres@vpWidthF = 0.60 ; Width of Labelbar
|
|
|
109 |
lres@vpHeightF = 0.10 ; Height of Labelbar
|
|
|
110 |
lres@lbPerimOn = False ; Turn off perimeter.
|
|
|
111 |
lres@lbOrientation = "Horizontal" ; Default is "Vertical"
|
|
|
112 |
lres@amLabelBarOrthogonalPosF = -1.0
|
|
|
113 |
lres@lbLabelAlignment = "InteriorEdges" ; Default is "BoxCenters".
|
|
|
114 |
lres@lbFillColors = cmap(2:,:) ; Colors for boxes.
|
|
|
115 |
lres@lbMonoFillPattern = True ; Fill them all solid.
|
|
|
116 |
lres@lbLabelFontHeightF = 0.015 ; Font Height
|
|
|
117 |
lres@lbLabelAutoStride = True ; Auto correct labels
|
|
|
118 |
lres@vpYF = 0.15 ; location of left edge
|
|
|
119 |
lres@vpXF = 0.22 ; location of top edge
|
|
|
120 |
|
|
|
121 |
|
|
|
122 |
; --- Draw polylines
|
|
|
123 |
do i = 0,ntra-1 ; Loop over all trajecories
|
|
|
124 |
|
|
|
125 |
x = traj(i,0:(ntim-1),1) ; Longitudes for trajectory i
|
|
|
126 |
y = traj(i,0:(ntim-1),2) ; Latitudes for trajectory i
|
|
|
127 |
|
|
|
128 |
do k = 0,ntim-2 ; Loop over all times
|
|
|
129 |
|
|
|
130 |
pres@gsLineColor = GetFillColor(cnLevels,cmap,traj(i,k,cnIndex)) ; color of polyline
|
|
|
131 |
gsn_polyline (wks,map,(/x(k),x(k+1)/),(/y(k),y(k+1)/),pres) ; draw polyline
|
|
|
132 |
|
|
|
133 |
end do
|
|
|
134 |
|
|
|
135 |
end do
|
|
|
136 |
|
|
|
137 |
|
|
|
138 |
; --- Draw colorbar/labelbar
|
|
|
139 |
gsn_labelbar_ndc(wks,dimsizes(cnLevels)+1,cnLabels,0.2,0.2,lres)
|
|
|
140 |
|
|
|
141 |
|
|
|
142 |
; ----------------- Markers ----------------------------------------------------------------
|
|
|
143 |
|
|
|
144 |
; --- Resources for the markers
|
|
|
145 |
kres = True ; marker resource
|
|
|
146 |
kres@gsMarkerSizeF = 5.0 ; marker size
|
|
|
147 |
kres@gsMarkerColor = "black" ; marker color
|
|
|
148 |
|
|
|
149 |
|
|
|
150 |
; --- Draw markers
|
|
|
151 |
do i = 0,ntra-1 ; Loop over all trajecories
|
|
|
152 |
|
|
|
153 |
x = traj(i,0:(ntim-1),1) ; Longitudes for trajectory i
|
|
|
154 |
y = traj(i,0:(ntim-1),2) ; Latitudes for trajectory i
|
|
|
155 |
|
|
|
156 |
do k = 0,ntim-2,5 ; Loop over all times
|
|
|
157 |
|
|
|
158 |
gsn_polymarker(wks,map,x(k),y(k),kres) ; mark position
|
|
|
159 |
|
|
|
160 |
end do
|
|
|
161 |
|
|
|
162 |
end do
|
|
|
163 |
|
|
|
164 |
|
|
|
165 |
; ----------------- Advance frame ----------------------------------- --------------------
|
|
|
166 |
|
|
|
167 |
frame(wks)
|
|
|
168 |
|