source: utils/tools_CPL_OASIS_ticket2379_2020/OBSTOOLS/dataplot/dataplot_txttimeseries.pro

Last change on this file was 3002, checked in by djlea, 9 years ago

Update documentation for obstools and dataplot. Removal of dataplot code not needed. Addition of headers to some dataplot code. Addition of .exe to command example in obstools.

File size: 10.9 KB
Line 
1pro plotts1, arrsv, title, typestr, minperc=minperc, $
2   juldatemin=juldatemin, juldatemax=juldatemax, $
3        emax=emax, emin=emin
4;+--------------------------------------------------------
5; plot mean and rms timeseries
6;
7; Author:  D. J. Lea      Feb 2008
8;+--------------------------------------------------------
9
10;date_label = LABEL_DATE(DATE_FORMAT = $
11;   ['%D %M, %Y'])
12
13;date_label = LABEL_DATE(DATE_FORMAT = $
14;   ['%D', '%M, %Y'])
15;date_label = LABEL_DATE(DATE_FORMAT = $
16;   ['%M %Y'])
17
18date_label = LABEL_DATE(DATE_FORMAT=['%D-%M','%Y'])
19
20
21; sort times (in case of a repeated day)
22
23timsrt=sort(arrsv(0,*))
24
25taxis=arrsv(0,timsrt)
26num=arrsv(1,timsrt)
27yaxis=arrsv(2,timsrt)
28yaxis2=arrsv(3,timsrt)
29
30; remove any zero times or non-finite values
31wh=where(taxis gt 0 and finite(yaxis) and finite(yaxis2))
32if (wh(0) gt -1) then begin
33taxis=taxis(wh)
34num=num(wh)
35yaxis=yaxis(wh)
36yaxis2=yaxis2(wh)
37endif
38
39; remove any with num lt than a specific value
40if (n_elements(minperc) eq 1) then begin
41maxnum=max(num,min=minnum)
42wh=where(num gt maxnum*minperc)
43if (wh(0) gt -1) then begin
44taxis=taxis(wh)
45num=num(wh)
46yaxis=yaxis(wh)
47yaxis2=yaxis2(wh)
48endif
49endif
50
51mxt=max(taxis,min=mnt)
52print, 'mnt mxt ',mnt, mxt
53
54ymx=max([yaxis,yaxis2],min=ymn)
55
56print, 'ymn ymx ',ymn,ymx
57
58;create a small amount of space around the max and min
59
60spc=(ymx-ymn)*0.05
61
62ymn=ymn-spc*3
63ymx=ymx+spc
64
65if (n_elements(emax) gt 0) then ymx=emax
66if (n_elements(emin) gt 0) then ymn=emin
67
68; setup time axis range
69
70skip=0
71xmx=max(taxis,min=xmn)
72if (n_elements(juldatemin) gt 0) then begin
73   if (xmn le juldatemin) then xmn=juldatemin
74        if (xmx le juldatemin) then skip=1
75endif
76if (n_elements(juldatemax) gt 0) then begin
77   if (xmx ge juldatemax) then xmx=juldatemax
78        if (xmn ge juldatemin) then skip=1
79endif
80
81if (skip eq 0) then begin
82;plot, arrsv(0,timsrt), arrsv(2,timsrt), xstyle=1, linestyle=1, $
83;   yrange=[ymn,ymx], $
84;   ytitle=typestr, title=title, $
85;   XTICKFORMAT = ['LABEL_DATE'], $
86;   XTICKUNITS = ['Day'], $
87;   XTICKINTERVAL = 4
88
89;plot, arrsv(0,timsrt), arrsv(3,timsrt), xstyle=1, /noerase, $
90;   yrange=[ymn,ymx], $
91;   XTICKFORMAT = ['LABEL_DATE'], $
92;   XTICKUNITS = ['Day'], $
93;   XTICKINTERVAL = 4
94
95;plot, arrsv(0,timsrt), arrsv(2,timsrt), xstyle=1, linestyle=1, $
96;   yrange=[ymn,ymx], $
97;   ytitle=typestr, title=title, $
98;   XTICKFORMAT = ['LABEL_DATE','LABEL_DATE'],$
99;   XTICKUNITS = ['Day','Month']
100
101;plot, arrsv(0,timsrt), arrsv(3,timsrt), xstyle=1, /noerase, $
102;   yrange=[ymn,ymx], $
103;   XTICKFORMAT = ['LABEL_DATE','LABEL_DATE'],$
104;   XTICKUNITS = ['Day','Month']
105
106;plot, taxis, yaxis, xstyle=1, linestyle=1, $
107;   yrange=[ymn,ymx], $
108;   ytitle=typestr, title=title, $
109;   XTICKFORMAT = ['LABEL_DATE']
110;   
111;plot, taxis, yaxis2, xstyle=1, /noerase, $
112;   yrange=[ymn,ymx], $
113;   XTICKFORMAT = ['LABEL_DATE']
114
115plot, taxis, yaxis, xstyle=1, ystyle=1, linestyle=1, $
116   yrange=[ymn,ymx], xrange=[xmn,xmx], $
117   ytitle=typestr, title=title, $
118   XTICKUNITS=['Time', 'Time'], XTICKFORMAT = ['LABEL_DATE'],YMARGIN=[6,4]
119
120plot, taxis, yaxis2, xstyle=4+1, ystyle=4+1, /noerase, $
121   yrange=[ymn,ymx], xrange=[xmn,xmx], $
122   XTICKUNITS=['Time', 'Time'], XTICKFORMAT = ['LABEL_DATE'],YMARGIN=[6,4]
123   
124
125
126; key
127      xcoord=0.8
128      ycoord=0.9
129      ycoord=0.35
130                ycoord=0.2
131
132                  ycoord2=ycoord-0.05
133                  xcoord2=xcoord+0.03
134                  xcoord3=xcoord+0.05
135                  plots, [xcoord,xcoord2],[ycoord,ycoord], linestyle=0, /normal
136                  xyouts, xcoord3, ycoord, 'RMS', /normal
137                  plots, [xcoord,xcoord2],[ycoord2,ycoord2], linestyle=1, /normal
138                  xyouts, xcoord3, ycoord2, 'mean',/normal               
139
140endif
141
142end
143
144
145; plot number
146
147pro plotts2, arrsv, title, typestr, minperc=minperc, $
148   juldatemin=juldatemin, juldatemax=juldatemax
149   
150
151; number
152;
153
154;date_label = LABEL_DATE(DATE_FORMAT = $
155;   ['%D %M, %Y'])
156
157;date_label = LABEL_DATE(DATE_FORMAT = $
158;   ['%M %Y'])
159
160date_label = LABEL_DATE(DATE_FORMAT=['%D-%M','%Y'])
161
162
163
164timsrt=sort(arrsv(0,*))
165
166taxis=arrsv(0,timsrt)
167yaxis=arrsv(1,timsrt)
168
169wh=where(taxis gt 0 and finite(yaxis))  ; remove any zero times and non-finite vals
170if (wh(0) gt -1) then begin
171taxis=taxis(wh)
172yaxis=yaxis(wh)
173endif
174mxt=max(taxis,min=mnt)
175print, 'mnt mxt ',mnt, mxt
176
177;plot, arrsv(0,timsrt), arrsv(1,timsrt), xstyle=1, $
178;   ytitle='Number of obs assim', title=title, $
179;   XTICKFORMAT = ['LABEL_DATE'], $
180;   XTICKUNITS = ['Day'], $
181;   XTICKINTERVAL = 4
182
183;info, taxis
184;info, yaxis
185
186;print,taxis, yaxis
187
188;plot, taxis, yaxis, xstyle=1, $
189;   ytitle='Number of obs assim', title=title, $
190;   XTICKFORMAT = ['LABEL_DATE']
191
192ymx=max(yaxis)*1.05
193
194; setup time axis range
195
196skip=0
197xmx=max(taxis,min=xmn)
198if (n_elements(juldatemin) gt 0) then begin
199   if (xmn le juldatemin) then xmn=juldatemin
200        if (xmx le juldatemin) then skip=1
201endif
202if (n_elements(juldatemax) gt 0) then begin
203   if (xmx ge juldatemax) then xmx=juldatemax
204        if (xmn ge juldatemin) then skip=1
205endif
206
207if (skip eq 0) then $
208plot, taxis, yaxis, xstyle=1, ystyle=1, $
209   ytitle='Number of obs assim', title=title, yrange=[0,ymx], xrange=[xmn,xmx],$
210   XTICKUNITS = ['Time', 'Time'], XTICKFORMAT = ['LABEL_DATE'],YMARGIN=[6,4]
211
212
213print,'min time ',min(arrsv(0,timsrt)),max(arrsv(0,timsrt))
214
215end
216
217PRO dataplot_txttimeseries, files, gif=gif, ps=ps, filtstr=filtstr, view=view, $
218   bin=bin, minperc=minperc, datemin=datemin, datemax=datemax, notitle=notitle, $
219        emax=emax, emin=emin
220
221; DJL switch off wave compatibility mode
222res=execute("waveoff")
223
224if (n_elements(filtstr) eq 0) then filtstr=""
225if (n_elements(view) eq 0) then view=0
226
227print, 'dataplot_txttimeseries '
228
229if (n_elements(datemin) gt 0) then begin
230   ; month, day, year
231   juldatemin=julday(datemin(1),datemin(2), datemin(0))
232        print, 'juldatemin set:', juldatemin, datemin
233endif
234if (n_elements(datemax) gt 0) then begin
235   ; month, day, year
236   juldatemax=julday(datemax(1),datemax(2), datemax(0))
237        print, 'juldatemax set:', juldatemax, datemax
238endif
239
240
241
242numfiles=n_elements(files)
243
244print,'numfiles ',numfiles
245
246imax=500
247   arrs=dblarr(4,imax)
248   arr=dblarr(4)
249   arrsv=dblarr(4,10000)
250
251j=0L     ; position in full array
252for ii=0,numfiles-1 do begin
253   print,files(ii)
254        OPENR, unit, files(ii), /get_lun
255        obstypestr=""
256        readf,unit,obstypestr
257        typestr=""
258        readf,unit,typestr
259        xrange=fltarr(2)
260        readf,unit,xrange
261        yrange=fltarr(2)
262        readf,unit,yrange
263        binspday=0.0
264        readf,unit,binspday
265
266   i=0
267   while (~ eof(unit) and i lt imax) do begin       
268        readf,unit,arr
269   print,i,arr
270        if arr(1) eq 0 then arr(2:*)=0
271        print,arr
272        arrs(*,i)=arr
273        i=i+1
274        endwhile
275 
276   numtimes=i
277   print, 'numtimes ',numtimes
278       
279        if (numtimes ge binspday) then begin
280         print,arrs(*,numtimes-binspday:numtimes-1)
281
282; store daily values from each file in full time series array
283
284   arrsv(*,j:j+binspday-1)=arrs(*,numtimes-binspday:numtimes-1)
285        j=j+binspday
286
287   endif else begin
288        if (numtimes gt 0) then begin
289       
290        print, '** numtimes ',numtimes
291        print, arrs(*,0:numtimes-1) 
292        arrsv(*,j:j+numtimes-1)=arrs(*,0:numtimes-1)       
293        j=j+numtimes
294               
295        endif
296        endelse
297               
298        FREE_LUN, unit
299       
300        print, obstypestr
301        print, typestr
302        print, xrange
303        print, yrange
304        print, binspday
305endfor   
306
307print, 'j ',j
308
309arrsv=arrsv(*,0:j-1)
310
311;stop
312
313; bin up the data
314
315;print,arrsv
316
317if (n_elements(bin) gt 0) then begin
318if (bin gt 1) then begin
319; time 0 num 1 mean 2 rms 3
320arrsv2=dblarr(4,j/bin)
321for i=0,j/bin-1 do begin
322
323   arrsvtemp=arrsv(*,i*bin:(i+1)*bin-1)
324        print,arrsvtemp
325   wh=where(arrsvtemp(1,*) gt 0)    ; number of obs gt 0
326;        print,wh
327        if (wh(0) gt -1) then begin
328
329; number of obs
330        arrsv2(1,i)=total(arrsvtemp(1,wh))
331; mean
332        arrsv2(2,i)=total(arrsvtemp(2,wh)*arrsvtemp(1,wh))/arrsv2(1,i)
333; rms
334        arrsv2(3,i)=sqrt(total(arrsvtemp(3,wh)^2*arrsvtemp(1,wh))/arrsv2(1,i))
335; date (average)
336;  print,arrsv(0,i*bin)
337   arrsv2(0,i)=total(arrsvtemp(0,*))/bin
338;  print,arrsv2(0,i)
339
340   endif
341
342endfor
343
344info, arrsv2
345
346arrsv=arrsv2
347   
348endif
349endif
350
351nel=n_elements(arrsv(0,*))
352
353; produce 1 month/3 month/all average values
354
355finaltime=arrsv(0,nel-1)
356onemon=finaltime-30
357threemon=finaltime-90
358
359print,arrsv
360
361wh1=where(arrsv(0,*) gt onemon and arrsv(1,*) gt 0)
362wh3=where(arrsv(0,*) gt threemon and arrsv(1,*) gt 0)
363wh0=where(arrsv(1,*) gt 0)
364
365num1=total(arrsv(1,wh1))
366mean1=total(arrsv(2,wh1)*arrsv(1,wh1))/num1
367rms1=sqrt(total(arrsv(3,wh1)^2*arrsv(1,wh1))/num1)
368print,num1
369
370num3=total(arrsv(1,wh3))
371mean3=total(arrsv(2,wh3)*arrsv(1,wh3))/num3
372rms3=sqrt(total(arrsv(3,wh3)^2*arrsv(1,wh3))/num3)
373print,num3
374
375num0=total(arrsv(1,wh0))
376mean0=total(arrsv(2,wh0)*arrsv(1,wh0))/num0
377rms0=sqrt(total(arrsv(3,wh0)^2*arrsv(1,wh0))/num0)
378print,num0
379
380if (keyword_set(notitle)) then begin
381title1=""
382title=""
383endif else begin
384subtitle=          'rms (mean), 1 month: '+strtrim(string(rms1,format='(G0.4)'),2)+$
385   '('+strtrim(string(mean1,format='(G0.3)'),2)+')'
386subtitle=subtitle+ ', 3 month: '+strtrim(string(rms3,format='(G0.4)'),2)+$
387   '('+strtrim(string(mean3,format='(G0.3)'),2)+')'
388subtitle=subtitle+ ', all: '+strtrim(string(rms0,format='(G0.4)'),2)+$
389   '('+strtrim(string(mean0,format='(G0.3)'),2)+')'
390
391
392fullfiltstr=''
393if (filtstr ne '') then fullfiltstr=' Type: '+filtstr
394
395title=obstypestr+typestr+'   lons ('+strtrim(string(xrange(0),format='(F0.2)'),2)+$
396   ','+strtrim(string(xrange(1),format='(F0.2)'),2)+$
397        ')   lats ('+strtrim(string(yrange(0),format='(F0.2)'),2)+','+$
398        strtrim(string(yrange(1),format='(F0.2)'),2)+')'+fullfiltstr
399
400title1=title+'!C'+subtitle
401endelse
402
403if (keyword_set(gif)) then begin
404   thisDevice = !D.Name
405
406        Set_Plot, 'Z'         ; do graphics in the background
407;  Device, Set_Resolution=[640,512], decomposed=0
408   Device, Set_Resolution=[800,512], decomposed=0
409        Erase                           ; clear any existing stuff
410        !p.charsize=0.75
411;        setupct, r, g, b     ; setup color table
412
413   plotts1, arrsv, title1, typestr, minperc=minperc, $
414   juldatemin=juldatemin, juldatemax=juldatemax, $
415        emax=emax, emin=emin
416
417
418   snapshot = TVRD()
419        WRITE_GIF,'dataplot_timeseries.gif',snapshot, r, g, b
420
421   plotts2, arrsv, title, typestr, minperc=minperc, $
422   juldatemin=juldatemin, juldatemax=juldatemax
423
424       
425   snapshot = TVRD()
426        WRITE_GIF,'dataplot_numtimeseries.gif',snapshot, r, g, b
427
428        Device, Z_Buffer=1    ; reset graphics mode
429        Set_Plot, thisDevice
430        !p.charsize=0.0
431
432endif
433
434if (keyword_set(ps)) then begin
435  ps=1
436  eps=0
437  landscape=1
438  pr2o,file='dataplot_timeseries.ps',landscape=landscape,ps=ps,eps=eps,color=1   
439  plotts1, arrsv, title1, typestr, minperc=minperc, $
440   juldatemin=juldatemin, juldatemax=juldatemax, $
441        emax=emax, emin=emin
442  prend2o,view=view
443 
444  pr2o,file='dataplot_numtimeseries.ps',landscape=landscape,ps=ps,eps=eps,color=1   
445  plotts2, arrsv, title, typestr, minperc=minperc, $
446   juldatemin=juldatemin, juldatemax=juldatemax
447  prend2o,view=view
448
449endif
450
451end
452
Note: See TracBrowser for help on using the repository browser.