New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
dataplot.pro in branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/TOOLS/OBSTOOLS/dataplot – NEMO

source: branches/2014/dev_r4650_UKMO14.5_SST_BIAS_CORRECTION/NEMOGCM/TOOLS/OBSTOOLS/dataplot/dataplot.pro @ 5967

Last change on this file since 5967 was 5967, checked in by timgraham, 8 years ago

Reset keywords before merging with head of trunk

  • Property svn:keywords set to Id
File size: 144.6 KB
Line 
1;+----------------------------------------------------------------------------------------
2; dataplot.pro
3; IDL widget based plotting routine for plotting observation and background values
4;
5; Author: D. J. Lea     -  Feb 2008
6;
7;+----------------------------------------------------------------------------------------
8;example calls:
9;
10;dataplot,filenamearray
11;dataplot,filenamearray,/batch,/gif       ; plots data directly to a gif
12;dataplot,filenamearray,/batch,/ps        ; plots data directly to a ps
13;dataplot,filenamearray,/batch,area=area            ; area 4 element array or descriptive string
14;dataplot, longitude, latitude, deparr, dayarr, valarr, bkgarr
15;dataplot, [longitude, latitude, deparr, dayarr, valarr, bkgarr], rmdi=rmdi, filename=filename
16;
17;optional keywords
18;area       descriptive string or array [minlon,minlat,maxlon,maxlat]
19;typeplot=1-5     plot obs-bkg or obs bkg etc
20;alldays=1     plot all days otherwise just plot first day   
21;depths        [depmin,depmax] in metres
22;showmdt    plot mdt
23;obstypeselect    string array of selected obs types to plot
24;printobstypes    keyword to print unique obs types to the terminal
25;
26;+-----------------------------------------------------------------------------------------
27; Note the upper two slider bars control the depth range
28; the lower two slider bars control the day range
29;------------------------------------------------------------------------------------------
30
31
32; Event-handling procedure.
33PRO dataplot_event, ev
34 
35  ; Retrieve the anonymous structure contained in the user value of
36  ; the top-level base widget. 
37 
38  WIDGET_CONTROL, ev.TOP, GET_UVALUE=stash
39  WIDGET_CONTROL, ev.ID, GET_UVALUE=uval
40
41;  print,uval
42
43;  print,'event'
44;  print,stash.depmin
45
46;  position = [!X.Window[0], !Y.Window[0], !X.Window[1], !Y.Window[1]]
47;  print,position
48;  xsize = (position(2) - position(0)) * !D.X_VSIZE
49;  ysize = (position(3) - position(1)) * !D.Y_VSIZE
50;  xstart = position(0) * !D.X_VSIZE
51;  ystart = position(1) * !D.Y_VSIZE
52;  print,xsize,ysize,xstart,ystart
53 
54  ; If the event is generated in the draw widget, update the
55  ; label values with the current cursor position and the value
56  ; of the data point under the cursor. Note that since we have
57  ; passed a pointer to the image array rather than the array
58  ; itself, we must dereference the pointer in the 'image' field
59  ; of the stash structure before getting the subscripted value.
60 
61  IF (TAG_NAMES(ev, /STRUCTURE_NAME) eq 'WIDGET_DRAW') THEN BEGIN
62   xdev= ev.X
63   ydev= ev.Y
64       
65;        print, 'xdev ',xdev, ' ydev ',ydev
66       
67; convert the pointer position to a data position
68   arr=CONVERT_COORD( xdev, ydev, /device, /to_data)
69   xdata=arr(0)
70   ydata=arr(1)
71
72;  print, xdata , ydata       
73       
74;        print, 'xdev: ',xdev,' ydev: ',ydev,' xdata: ',xdata, ' ydata: ',ydata
75;    WIDGET_CONTROL, stash.label1, $
76;      SET_VALUE='X position: ' + STRING(xdata)
77;    WIDGET_CONTROL, stash.label2, $
78;      SET_VALUE='Y position: ' + STRING(ydata)
79;;    WIDGET_CONTROL, stash.label3, $
80;;      SET_VALUE='Value: ' + $
81;;      STRING((*stash.imagePtr)[ev.X, ev.Y], FORMAT='(Z12)')
82;    WIDGET_CONTROL, stash.label3, $
83
84
85   ; What kind of event is this?
86
87;eventTypes = ['DOWN', 'UP', 'MOTION']
88;0 1 2
89;thisEvent = ev.type
90
91;print, 'ev.type ',ev.type
92
93CASE ev.type OF
94   0: BEGIN
95; dragbox start
96   print, 'down ',xdev,ydev, xdata, ydata
97         ; Turn motion events on for the draw widget.
98
99      Widget_Control, stash.draw, Draw_Motion_Events=1
100
101   ; dragbox
102         ; Create a pixmap. Store its ID. Copy window contents into it.
103
104      Window, /Free, /Pixmap, XSize=stash.im_size(1), YSize=stash.im_size(2)
105      stash.pixID = !D.Window
106      Device, Copy=[0, 0, stash.im_size(1), stash.im_size(2), 0, 0, stash.drawID]
107;      WSet, stash.drawID
108
109   print,stash.im_size(0), stash.im_size(1), stash.im_size(2), stash.pixID, stash.drawID
110     
111      stash.xcorn=xdev
112      stash.ycorn=ydev
113     
114      stash.xdatacorn=xdata
115      stash.ydatacorn=ydata
116   
117;   RETURN
118   ENDCASE
119   1: BEGIN
120; dragbox close
121   print, 'up ',xdev,ydev, xdata, ydata
122         ; Turn draw motion events off. Clear any events queued for widget.
123
124      Widget_Control, stash.draw, Draw_Motion_Events=0, Clear_Events=1
125
126       WSet, stash.drawID
127      Device, Copy=[0, 0, stash.im_size(1), stash.im_size(2), 0, 0, stash.pixID]
128      WDelete, stash.pixID
129           
130      ; zoom in
131
132; dragbox   
133; avoid zooming if no dragging has occurred
134   pixdiffx = min([abs(stash.xcorn-xdev),abs(stash.ycorn-ydev)])       
135        print, '** Pixdiffx ',pixdiffx, stash.xcorn-xdev, stash.ycorn-ydev
136        if (pixdiffx gt 1 and finite(xdata) and finite(ydata)) then begin
137
138      minxdata = min([stash.xdatacorn, xdata], max=maxxdata)
139        minydata = min([stash.ydatacorn, ydata], max=maxydata)
140
141
142   if (minxdata ge -180 and maxxdata le 360 and minydata ge -90 and maxydata le 90) then begin
143
144                stash.xrange=[minxdata,maxxdata]
145         stash.yrange=[minydata,maxydata]
146               
147           plotpoints,stash
148        endif
149        endif 
150           
151;   RETURN
152   ENDCASE
153   2: BEGIN
154; drag box motion
155;   print, 'motion: ',ev.x, ev.y
156
157      WSet, stash.drawID
158      Device, Copy=[0, 0, stash.im_size(1), stash.im_size(2), 0, 0, stash.pixID]
159
160   
161   sx = stash.xcorn
162   sy = stash.ycorn
163
164      PlotS, [sx, sx, xdev, xdev, sx], [sy, ydev, ydev, sy, sy], /Device, $
165          Color=!p.color
166;         Color=255
167   
168   ENDCASE
169   
170;   RETURN
171   7: begin
172      print,'mouse wheel event ',ev.clicks
173      if (ev.clicks eq 1 or ev.clicks eq -1) then begin
174      print,'zoom in/out'
175
176      print,'stash.depmin ',stash.depmin
177
178      print,'xdata ',xdata
179               
180                if (finite(xdata) and finite(ydata)) then begin
181
182; make an adjustment if we're in the Pacific
183   xrange=stash.xrange
184        typeproj=stash.typeproj
185   if (xrange(1) gt 180 and typeproj eq 1) then if (xdata lt 0) then xdata=xdata+360
186       
187
188      oxrange=stash.xrange
189      oyrange=stash.yrange
190
191      dataplot_zoom, stash, xdata, ydata, ev.clicks
192
193; plot only if a change has been made
194      if (stash.xrange(0) ne oxrange(0) or stash.yrange(0) ne oyrange(0) $
195                    or stash.xrange(1) ne oxrange(1) or stash.yrange(1) ne oyrange(1)) then $
196         plotpoints, stash
197
198      endif
199               
200                endif
201
202   ENDCASE
203ENDCASE
204
205
206
207
208
209
210
211
212;      SET_VALUE='Value: ' + STRING(ev.release)
213
214;  if(ev.release eq 1) then begin
215;     print,'left click - search for nearest point'
216   if(ev.release eq 1 or ev.release eq 4) then begin     ; mouse click
217
218; might want a generalised point selection routine
219; for this routine and for plotpoints
220
221      print,'xdata ',xdata, ' ydata ',ydata
222                print,'xdev ',xdev, 'ydev ',ydev
223
224      if (finite(xdata) and finite(ydata)) then begin
225
226      print, 'finite'
227
228      xarr1=stash.xarr
229      yarr1=stash.yarr
230      dep1=stash.dep
231      dayarr=stash.dayarr
232      daymin=stash.daymin
233      daymax=stash.daymax
234                obstypes1=stash.obstypes
235      xrange=stash.xrange
236      yrange=stash.yrange
237      xs=xrange(1)-xrange(0)
238      ys=yrange(1)-yrange(0)
239      rmdi=stash.rmdi
240;                if (stash.salinity eq 0) then begin
241;       obs1=stash.obs
242;       bkg1=stash.bkg
243;                 qcarr1=stash.qcarr
244;                endif else begin
245;       obs1=stash.obs2
246;       bkg1=stash.bkg2
247;                 qcarr1=stash.qcarr2
248;     endelse                   
249                typeproj=stash.typeproj
250
251; make an adjustment if we're in the Pacific
252   if (xrange(1) gt 180 and typeproj eq 1) then if (xdata lt 0) then xdata=xdata+360
253
254; set reasonable dist criteria
255
256      distmax=(xs/100.)^2+(ys/100.)^2
257;     print,'distmax ',distmax
258
259
260      selpoints,stash,lonsel, latsel, innovsel, qcsel, daysel, obstypsel, obnumsel, numsel, typestr
261
262; get max min daysel
263      mindaysel=min(daysel, max=maxdaysel)
264
265;     print, 'selpoints ',obnumsel
266
267      if (n_elements(lonsel) gt 0) then begin
268      dist=((lonsel-xdata)^2+(latsel-ydata)^2)
269      res=min(dist,whd)    ; whd is the minumum dist point index
270
271;     print,lonsel
272
273   xselstr=""
274   yselstr=""
275   valselstr=""
276   qcselstr=""
277   datestr=""
278        obnumselstr=""
279        mindatestr=""
280        maxdatestr=""
281
282;  print,'whd(0) ',whd(0)
283;        print,'res ',res, ' distmax ',distmax
284
285   obstypstr=""     
286   if (whd(0) gt -1 and res lt distmax) then begin   
287   xsel=lonsel(whd)
288   ysel=latsel(whd)
289 
290   print,'xsel ',xsel,' ysel ',ysel       
291 ; make an adjustment if we're in the Pacific
292   if (xrange(1) gt 180 and typeproj eq 1) then if (xsel gt 180) then xsel=xsel-360
293               
294   valsel=innovsel(whd)
295   qcsel=qcsel(whd)
296   datesel=daysel(whd)
297   obstypsel=obstypsel(whd)
298   obnumsel=obnumsel(whd)
299
300   xselstr=string(xsel)
301   yselstr=string(ysel)
302   valselstr=string(valsel)   
303   qcselstr=string(qcsel)
304;  datestr=string(datesel)
305   obstypstr=string(obstypsel(0))
306        obnumselstr=string(obnumsel(0))
307        print,'datesel(0) ',datesel(0)
308        print,'obnumsel(0) ',obnumsel(0)
309   jul_to_dtstr,datesel(0),datestr
310        jul_to_dtstr,mindaysel,mindatestr
311        jul_to_dtstr,maxdaysel,maxdatestr
312   
313   if (ev.release eq 4) then begin
314       
315         print,'ev.release ',ev.release
316
317;     window
318      wh=where(xarr1 eq xsel and yarr1 eq ysel)
319                print,'xsel ',xsel
320;                print,'xarr1 ',xarr1
321      dep2=dep1(wh)
322
323      if (wh(0) gt -1) then begin
324      if (stash.density eq 0) then begin
325
326         if (stash.salinity eq 0) then begin
327;                  print,'salinity eq 0'
328            obs2=stash.obs(wh)
329            bkg2=stash.bkg(wh)
330            qcarr2=stash.qcarr(wh)
331         endif else begin
332;                  print,'salinity eq 1'
333            obs2=stash.obs2(wh)
334            bkg2=stash.bkg2(wh)
335            qcarr2=stash.qcarr2(wh)
336         endelse
337
338      endif else begin
339
340         if (stash.filetype eq 'CRT') then begin
341
342            bkgU=stash.bkg(wh)
343            bkgV=stash.bkg2(wh)
344            bkg2=sqrt(bkgU*bkgU + bkgV*bkgV)
345
346            obsU=stash.obs(wh)
347            obsV=stash.obs2(wh)
348            obs2=sqrt(obsU*obsU + obsV*obsV)
349
350            qcarr2=stash.qcarr(wh)
351
352         endif else begin
353
354;            print,'density eq 1'
355            bkgt=stash.bkg(wh)
356            bkgs=stash.bkg2(wh)
357            obst=stash.obs(wh)
358            obss=stash.obs2(wh)                 
359;            obs2=obst+obss
360;            bkg2=bkgt+bkgs
361            obs2=eos_nemo(obst,obss)
362            bkg2=eos_nemo(bkgt,bkgs)
363            wh2=where(bkgt eq rmdi or bkgs eq rmdi)
364            if (wh2(0) gt -1) then begin
365               obs2(wh2)=rmdi
366               bkg2(wh2)=rmdi
367            endif
368            qcarr2=max([[stash.qcarr(wh)],[stash.qcarr2(wh)]],dim=2)                                     
369
370         endelse
371
372      endelse
373
374;     val2=abs(obs1(wh)-bkg1(wh))
375      val2=abs(obs2-bkg2)
376;     obs2=obs1(wh)
377;     bkg2=bkg1(wh)
378      obstype2=obstypes1(wh)
379;                qcarr2=qcarr1(wh)
380                dayarr2=dayarr(wh)
381;     val2=sqrt(val1(wh)^2)
382     
383;     print,'val2 ',val2(wh2)
384;     print,'dep2 ',dep2(wh2)
385;     plot,val2(wh2),dep2(wh2),ystyle=1,xstyle=1
386
387      profilewindow,dep2,val2,obs2,bkg2,obstype2,qcarr2,dayarr2,datestr,xselstr,yselstr,rmdi, salinity=stash.salinity, $
388                  plot_bad_obs=stash.plot_bad_obs, white_on_black=stash.white_on_black, coltable=stash.coltable, $
389                           depmax=stash.depmax, depmin=stash.depmin
390
391
392; set graphics back to main window
393;        WSET, stash.drawID
394      print,'ev release'
395      plotpoints,stash
396               
397                endif     ; wh(0)>-1
398
399   endif       ; ev.release 4 (right mouse)
400        endif         
401   endif       ; left or right mouse click
402   
403    WIDGET_CONTROL, stash.label1, $
404      SET_VALUE=stash.labt1 + strtrim(xselstr,1)
405    WIDGET_CONTROL, stash.label2, $
406      SET_VALUE=stash.labt2 + strtrim(yselstr,1)
407     WIDGET_CONTROL, stash.label3, $
408      SET_VALUE=typestr+': ' + strtrim(valselstr,1)
409     WIDGET_CONTROL, stash.label4, $
410      SET_VALUE=stash.labt4 + strtrim(qcselstr,1)
411     WIDGET_CONTROL, stash.label5, $
412      SET_VALUE=stash.labt5 + strtrim(datestr,1)
413     WIDGET_CONTROL, stash.label6, $
414      SET_VALUE=stash.labt6 + strtrim(valselstr,1)
415;     WIDGET_CONTROL, stash.label7, $
416;      SET_VALUE=stash.labt7 + strtrim(string(innovsd),1)
417     WIDGET_CONTROL, stash.label8, $
418      SET_VALUE=stash.labt8 + strtrim(obstypstr,1)
419     WIDGET_CONTROL, stash.label9, $
420      SET_VALUE=stash.labt9 + strtrim(obnumselstr,1)
421     WIDGET_CONTROL, stash.label10, $
422      SET_VALUE=stash.labt10 + string(mindatestr) + ' ' + string(maxdatestr)
423     
424   endif else begin        ; not finite
425
426   print, 'not finite'
427       
428; check for colorbar click
429      if (ydev lt 80) then begin       
430           toplevel=stash.base
431           fmx=stash.fmx
432         fmn=stash.fmn
433                 mx=stash.mx
434         mn=stash.mn
435         inputmxmnwindow, fmx, fmn, mx, mn, stash.rmdi, toplevel, success
436      if (success eq 1) then begin
437         stash.fmx=fmx
438         stash.fmn=fmn
439         plotpoints,stash
440         endif
441                endif
442       
443        endelse         ; finite
444
445   endif
446
447     
448     
449;      print,tag_names(ev)
450;      print,ev.id
451;      print,ev.top
452;      print,ev.handler
453;      print,ev.type
454;      print,ev.press
455;      print,ev.release
456;      print,ev.clicks
457;      print,ev.modifiers
458;      print,ev.ch
459;      print,ev.key
460       
461   ENDIF                   ; end of widget draw events
462 
463  ; If the event is generated in a button, destroy the widget
464  ; hierarchy. We know we can use this simple test because there
465  ; is only one button in the application.
466 
467
468;  IF (TAG_NAMES(ev, /STRUCTURE_NAME) eq 'WIDGET_COMBOBOX') THEN BEGIN
469;   IF (uval eq 'LEVELLIST') THEN BEGIN
470
471   IF (uval eq 'LEVELCHOICE') THEN BEGIN
472
473   stash.depmin=ev.value
474
475      print,'levelchoice'
476   print, 'stash.depmin: ',stash.depmin
477   print, 'stash.depmax: ',stash.depmax
478
479   if (stash.depmin ge stash.depmax) then begin
480      stash.depmax=stash.depmin
481   if (stash.depmax gt stash.depmaxl) then stash.depmax=stash.depmaxl
482
483;set sliderb position
484   WIDGET_CONTROL, stash.sliderb, set_value=stash.depmax
485
486   endif
487   
488
489      plotpoints, stash
490
491     
492;        print,ev.index,widget_info(ev.id,/combobox_gettext)
493 
494  ENDIF
495
496   IF (uval eq 'LEVELCHOICEB') THEN BEGIN
497
498   stash.depmax=ev.value
499
500      print,'levelchoiceb'
501   print, 'stash.depmin: ',stash.depmin
502   print, 'stash.depmax: ',stash.depmax
503
504   if (stash.depmax le stash.depmin) then begin
505      stash.depmin=stash.depmax
506   if (stash.depmin lt stash.depminl) then stash.depmin=stash.depminl
507
508;set sliderb position
509   WIDGET_CONTROL, stash.slider, set_value=stash.depmin
510
511   endif
512   
513
514      plotpoints, stash
515
516     
517;        print,ev.index,widget_info(ev.id,/combobox_gettext)
518 
519  ENDIF
520
521
522  IF (uval eq 'DATERANGE1') THEN BEGIN
523
524   odaymin=stash.daymin
525        odaymax=stash.daymax
526 
527   stash.daymin=ev.value
528;  print, 'ev.drag ',ev.drag
529
530; set slider label
531   jul_to_dtstr,stash.daymin,dayminstr, /notime
532        WIDGET_CONTROL, stash.slider1label, $
533        SET_VALUE='min date: '+dayminstr
534
535   if (stash.daymin ge stash.daymax) then begin
536;     stash.daymax=stash.daymin+1
537      stash.daymax=stash.daymin
538   if (stash.daymax gt stash.daymaxl) then stash.daymax=stash.daymaxl
539
540;set slider2 position
541   WIDGET_CONTROL, stash.slider2, set_value=stash.daymax
542; set slider 2 label
543   jul_to_dtstr,stash.daymax,daymaxstr, /notime
544        WIDGET_CONTROL, stash.slider2label, $
545        SET_VALUE='max date: '+daymaxstr
546   
547       
548   endif
549   
550   print,'daterange1'
551;  if (ev.drag eq 0) then
552        if (stash.daymin ne odaymin or stash.daymax ne odaymax) then plotpoints,stash
553   
554  ENDIF
555  IF (uval eq 'DATERANGE2') THEN BEGIN
556
557   odaymin=stash.daymin
558        odaymax=stash.daymax
559 
560   stash.daymax=ev.value
561
562; set slider label
563   jul_to_dtstr,stash.daymax,daymaxstr, /notime
564        WIDGET_CONTROL, stash.slider2label, $
565        SET_VALUE='max date: '+daymaxstr
566
567   if (stash.daymax le stash.daymin) then begin
568;     stash.daymin=stash.daymax-1
569      stash.daymin=stash.daymax
570   if (stash.daymin lt stash.dayminl) then stash.daymin=stash.dayminl
571
572;set slider1 position
573   WIDGET_CONTROL, stash.slider1, set_value=stash.daymin
574; set slider 2 label
575   jul_to_dtstr,stash.daymin,dayminstr, /notime
576        WIDGET_CONTROL, stash.slider1label, $
577        SET_VALUE='min date: '+dayminstr
578
579   endif
580   
581   print,'daterange2'
582;  if (ev.drag eq 0) then
583        if (stash.daymin ne odaymin or stash.daymax ne odaymax) then plotpoints,stash
584
585  ENDIF
586 
587   IF (uval eq "RADIO1") THEN BEGIN
588      print,'uval eq RADIO1'
589      stash.typeplot=1
590                print, 'typeplot ',stash.typeplot
591      plotpoints,stash
592   ENDIF
593   IF (uval eq "RADIO2") THEN BEGIN
594      print,'uval eq RADIO2' 
595      stash.typeplot=2
596                print, 'typeplot ',stash.typeplot
597      plotpoints,stash
598   ENDIF
599   IF (uval eq "RADIO3") THEN BEGIN
600      print,'uval eq RADIO3'
601      stash.typeplot=3
602                print, 'typeplot ',stash.typeplot
603      plotpoints,stash
604   ENDIF
605   IF (uval eq "RADIO4") THEN BEGIN
606      print,'uval eq RADIO4'
607      stash.typeplot=4
608                print, 'typeplot ',stash.typeplot
609      plotpoints,stash
610   ENDIF
611
612   IF (uval eq "RADIO5") THEN BEGIN
613      print,'uval eq RADIO5'
614      stash.ombtypeplot=1
615                print, 'ombtypeplot ',stash.ombtypeplot
616      plotpoints,stash
617   ENDIF
618   IF (uval eq "RADIO6") THEN BEGIN
619      print,'uval eq RADIO6'
620      stash.ombtypeplot=2
621                print, 'ombtypeplot ',stash.ombtypeplot
622      plotpoints,stash
623   ENDIF
624   IF (uval eq "RADIO7") THEN BEGIN
625      print,'uval eq RADIO7'
626      stash.ombtypeplot=3
627                print, 'ombtypeplot ',stash.ombtypeplot
628      plotpoints,stash
629   ENDIF
630
631   IF (uval eq "RADIO01") THEN BEGIN
632      print,'uval eq RADIO01'
633      stash.typeproj=1
634      stash.xrange=stash.xrangedef
635      stash.yrange=stash.yrangedef
636      plotpoints,stash
637   ENDIF
638   IF (uval eq "RADIO02") THEN BEGIN
639      print,'uval eq RADIO02'
640      stash.typeproj=2
641      stash.xrange=stash.xrangedef
642      stash.yrange=stash.yrangedef
643      plotpoints,stash
644   ENDIF
645   IF (uval eq "RADIO03") THEN BEGIN
646      print,'uval eq RADIO03'
647      stash.typeproj=3
648      stash.xrange=stash.xrangedef
649      stash.yrange=stash.yrangedef
650      plotpoints,stash
651   ENDIF
652
653   IF (uval eq "RADIO001") THEN BEGIN
654      print,'uval eq RADIO001 ev.select: ',ev.select
655      stash.plot_bad_obs=ev.select
656;     stash.xrange=stash.xrangedef
657;     stash.yrange=stash.yrangedef
658      plotpoints,stash
659   ENDIF
660
661
662   IF (uval eq "RADIOSAL") THEN BEGIN
663       print,'uval eq RADIOSAL ev.select: ',ev.select
664       salinity=ev.select
665;            if (stash.filetype eq "Prof" or stash.filetype eq "feedback") then begin
666          stash.salinity=salinity
667          plotpoints,stash
668;      endif
669   ENDIF
670
671
672   IF (uval eq "RADIODENSITY") THEN BEGIN
673            print,'uval eq RADIODENSITY ev.select: ',ev.select
674            density=ev.select
675            if (stash.filetype eq "CRT") then stash.salinity=0
676;            if (stash.filetype eq "Prof" or stash.filetype eq "CRT") then begin
677               stash.density=density 
678               plotpoints,stash
679;           endif
680   ENDIF
681
682
683  IF (uval eq "PRINT") THEN BEGIN
684   ps=1
685   eps=0
686   landscape=1
687   pr2,file=stash.outfile+'.ps',landscape=landscape,ps=ps,eps=eps,color=1 
688   plotpoints,stash
689   prend2,/view
690  ENDIF
691
692  IF (uval eq "SAVE") THEN BEGIN
693
694   thisDevice = !D.Name
695        psave=!p
696
697        Set_Plot, 'Z'         ; do graphics in the background
698;  Device, Set_Resolution=[800,512], decomposed=0
699        Device, Set_Resolution=stash.picsize, decomposed=0
700        Erase                           ; clear any existing stuff
701        !p.charsize=0.75
702
703;  if (stash.white_on_black eq 0) then begin
704;; flip background and foreground color
705;                        pcolor=!p.color
706;                        pbackground=!p.background
707;                        !p.color=pbackground
708;                        !p.background=pcolor
709;
710;        endif
711;        print,'!p.color,!p.background ',!p.color,!p.background
712
713        setupct, r, g, b, coltable=stash.coltable, $
714         white_on_black=stash.white_on_black    ; setup color table
715; plot data
716   plotpoints,stash
717   snapshot = TVRD()
718        WRITE_GIF,stash.outfile+'.gif',snapshot, r, g, b
719        Device, Z_Buffer=1    ; reset graphics mode
720        Set_Plot, thisDevice
721;        !p.charsize=0.0
722   !p=psave
723       
724        spawn,'xv '+stash.outfile+'.gif'
725
726  ENDIF
727
728;Area selection
729
730  xrange=stash.xrange
731  yrange=stash.yrange
732  areasel,uval,xrange,yrange,success, toplevel=stash.base
733  if (success eq 1) then begin
734  stash.xrange=xrange
735  stash.yrange=yrange
736  plotpoints,stash
737  endif
738
739     
740  IF (uval eq "Timeseries") THEN BEGIN
741   stash.numtimeseries=0
742        timeserieswindow,stash
743; set graphics back to main window
744;        WSET, stash.drawID
745;;       print,'ev release'
746;;    plotpoints,stash
747
748  ENDIF
749
750  IF (uval eq "Num timeseries") THEN BEGIN
751   stash.numtimeseries=1 
752        timeserieswindow,stash
753; set graphics back to main window
754;        WSET, stash.drawID
755      print,'ev release'
756      plotpoints,stash
757
758  ENDIF
759
760  IF (uval eq "TS diagram") THEN BEGIN
761        tsdiagramwindow,stash
762; set graphics back to main window
763;     WSET, stash.drawID
764   print,'ev release'
765   plotpoints,stash
766
767  ENDIF
768
769  IF (uval eq "Worst points") THEN BEGIN
770        worstpointswindow,stash
771  ENDIF
772
773  IF (uval eq "Filter type") THEN BEGIN
774   filterwindow, stash
775  ENDIF
776
777  IF (uval eq "Low res map") THEN BEGIN
778;        stash.map_file=""
779        stash.hires_map=0       
780        Widget_Control, stash.loresmap, Set_Button=1
781        Widget_Control, stash.hiresmap, Set_Button=0
782        plotpoints,stash
783  ENDIF
784
785  IF (uval eq "Hi res map") THEN BEGIN
786;        stash.map_file="/opt/ukmo/wave/ukmo/data/map_europe.xdr"
787   stash.hires_map=1
788        Widget_Control, stash.loresmap, Set_Button=0
789        Widget_Control, stash.hiresmap, Set_Button=1
790        plotpoints,stash
791  ENDIF
792
793  IF (uval eq "Plot only bad obs") THEN BEGIN
794   stash.plot_only_bad_obs = 1-stash.plot_only_bad_obs
795        Widget_Control, stash.pltonbado, Set_Button=stash.plot_only_bad_obs
796        if (stash.plot_bad_obs eq 1) then plotpoints,stash
797  ENDIF
798
799  IF (uval eq "White on black") THEN BEGIN
800        print,uval,' !p.color,!p.background ',!p.color,!p.background
801   stash.white_on_black = 1-stash.white_on_black
802        Widget_Control, stash.whtonblack, Set_Button=stash.white_on_black
803;        IF (!d.name eq 'X') then begin
804;                 WSET, stash.drawID
805;; flip background and foreground color
806;                        pcolor=!p.color
807;                        pbackground=!p.background
808;                        !p.color=pbackground
809;                        !p.background=pcolor
810
811        setupct, r, g, b, coltable=stash.coltable, $
812         white_on_black=stash.white_on_black    ; setup color table
813         
814                        plotpoints,stash
815;  endif
816
817
818  ENDIF
819
820
821  IF (uval eq "Square psym") THEN BEGIN
822        stash.sym=1
823        plotpoints,stash
824  ENDIF
825  IF (uval eq "Round psym") THEN BEGIN
826        stash.sym=2
827        plotpoints,stash
828  ENDIF
829
830  IF (uval eq "incsym") THEN BEGIN
831   stash.symscale=stash.symscale*1.41421
832        plotpoints,stash
833  ENDIF
834  IF (uval eq "decsym") THEN BEGIN
835   stash.symscale=stash.symscale/1.41421
836        plotpoints,stash
837  ENDIF
838  IF (uval eq "resetsym") THEN BEGIN
839   stash.symscale=1.0
840        plotpoints,stash
841  ENDIF
842
843  IF (uval eq "vertgrad") THEN BEGIN
844   stash.vertgrad=1-stash.vertgrad
845        Widget_Control, stash.vertgradmenu, Set_Button=stash.vertgrad
846        plotpoints,stash
847  ENDIF
848
849  IF (uval eq "Info") THEN BEGIN
850         infowindow
851  ENDIF
852
853  IF (uval eq "input max/min") THEN BEGIN
854        toplevel=stash.base
855        fmx=stash.fmx
856        fmn=stash.fmn
857        mx=stash.mx
858        mn=stash.mn
859      inputmxmnwindow, fmx, fmn, mx, mn, stash.rmdi, toplevel, success
860   if (success eq 1) then begin
861        stash.fmx=fmx
862        stash.fmn=fmn
863        plotpoints,stash
864        endif
865       
866;        stash.mx=mx
867;        stash.mn=mn     
868;        info,stash.fmx
869;        info,stash.fmn
870;        info,stash.mx
871;        info,stash.mn
872
873  ENDIF
874
875     
876; store values
877  WIDGET_CONTROL, ev.TOP, SET_UVALUE=stash
878
879;  IF (TAG_NAMES(ev, /STRUCTURE_NAME) eq 'WIDGET_BUTTON') THEN BEGIN
880  IF (uval eq "DONE") THEN BEGIN
881;    WIDGET_CONTROL, ev.TOP, /DESTROY
882    WIDGET_CONTROL, stash.base, /DESTROY
883  ENDIF
884
885
886
887;  WIDGET_CONTROL, stash.label3, SET_VALUE=TAG_NAMES(ev, /STRUCTURE_NAME)
888;  print,TAG_NAMES(ev, /STRUCTURE_NAME), uval
889 
890END
891
892;-----------------
893; proceedures
894;-----------------
895
896;setup area ranges
897
898PRO areasel, uval, xrange, yrange, success, toplevel=toplevel
899
900;info, xrange
901;info, yrange
902
903success=0
904 
905  IF (uval eq "Global") THEN BEGIN
906      xrange=[-180.,180.]
907        yrange=[-90.,90.]
908        success=1
909  ENDIF
910  IF (uval eq "Arctic") THEN BEGIN
911      xrange=[-180.,180.]
912        yrange=[70.,90.]
913        success=1
914  ENDIF
915  IF (uval eq "N Atl") THEN BEGIN
916      xrange=[-100.,0.]
917        yrange=[25.,70.]
918        success=1
919  ENDIF
920  IF (uval eq "Trop Atl") THEN BEGIN
921      xrange=[-80.,20.]
922        yrange=[-25.,25.]
923        success=1
924  ENDIF
925  IF (uval eq "S Atl") THEN BEGIN
926      xrange=[-70.,20.]
927        yrange=[-50.,-25.]
928        success=1
929  ENDIF
930  IF (uval eq "N Pac") THEN BEGIN     
931        xrange=[120.,-100.+360.]
932        yrange=[25.,70.]
933        success=1
934  ENDIF
935  IF (uval eq "Trop Pac") THEN BEGIN
936      xrange=[120.,-80.+360.]
937        yrange=[-25.,25.]
938        success=1
939  ENDIF
940  IF (uval eq "S Pac") THEN BEGIN
941      xrange=[120.,-70.+360.]
942        yrange=[-50.,-25.]
943        success=1
944  ENDIF
945  IF (uval eq "Indian") THEN BEGIN
946      xrange=[20.,120.]
947        yrange=[-50.,30.]
948        success=1
949  ENDIF
950  IF (uval eq "S Ocean") THEN BEGIN
951      xrange=[-180.,180.]
952        yrange=[-90.,-50.]
953        success=1
954  ENDIF
955  IF (uval eq "Pacific") THEN BEGIN
956      xrange=[120.,-80.+360.]
957        yrange=[-50.,70.]
958        success=1
959  ENDIF
960   IF (uval eq "Atlantic") THEN BEGIN
961      xrange=[-100.,20.]
962        yrange=[-50.,70.]
963        success=1
964  ENDIF
965   IF (uval eq "Med") THEN BEGIN
966      xrange=[-15.,38.]
967        yrange=[30.,46.]
968        success=1
969  ENDIF
970   IF (uval eq "NWS") THEN BEGIN
971      xrange=[-20.,13.]
972        yrange=[40.,65.]
973        success=1
974  ENDIF
975   IF (uval eq "Input Area") THEN BEGIN
976      IF (n_elements(toplevel) gt 0) THEN BEGIN
977        print,'success b4',success
978      inputareawindow, xrange, yrange, toplevel, success
979        print,'success af',success
980        ENDIF
981   ENDIF
982
983END
984
985 
986;select points to plot
987
988PRO selpoints, stash, lonsel, latsel, innovsel, qcsel, daysel, obstypsel, obnumsel, numsel, typestr, $
989   daymin=daymin, daymax=daymax, salinity=salinity, rmsmean=rmsmean, innovsel2=innovsel2
990
991;     print, 'calling selpoints ',stash.netcdf, stash.txt
992      print, 'calling selpoints'
993
994      xarr1=stash.xarr
995      yarr1=stash.yarr
996      dep1=stash.dep
997                if (n_elements(salinity) eq 0) then salinity=stash.salinity
998      density=stash.density
999                mld=stash.mld
1000      filetype=stash.filetype
1001
1002      obnum=stash.obnum
1003;     obnum=lindgen(n_elements(xarr1)) ; generate an array of observation numbers
1004;                             ; starting with zero
1005
1006;          if (salinity eq 0) then begin                 
1007;       obs1=stash.obs
1008;       bkg1=stash.bkg
1009;       qcarr=stash.qcarr
1010;                endif else begin
1011;          obs1=stash.obs2
1012;       bkg1=stash.bkg2
1013;       qcarr=stash.qcarr2
1014;                endelse             
1015
1016      if (salinity eq 0) then qcarr=stash.qcarr
1017      if (salinity eq 1 or mld eq 1) then qcarr=stash.qcarr2
1018      if (density  eq 1) then begin
1019         if (filetype eq 'CRT') then qcarr=stash.qcarr $
1020         else qcarr=max([[stash.qcarr],[stash.qcarr2]],dim=2)
1021      endif
1022
1023      obstypes1=stash.obstypes
1024      depmin=stash.depmin
1025      depmax=stash.depmax
1026      xrange=stash.xrange
1027      yrange=stash.yrange
1028      dayarr=stash.dayarr
1029                fdayarr=long(dayarr+stash.dayshi)  ; 0:00 hours should be plotted in prev day
1030      if (n_elements(daymin) eq 0) then daymin=stash.daymin
1031      if (n_elements(daymax) eq 0) then daymax=stash.daymax
1032      fdaymin=long(daymin)
1033                fdaymax=long(daymax)
1034      typeplot=stash.typeplot
1035                ombtypeplot=stash.ombtypeplot
1036      typeproj=stash.typeproj
1037      xrange=stash.xrange
1038      yrange=stash.yrange
1039      rmdi=stash.rmdi
1040
1041      print,'dayarr ',max(dayarr)-2454710d, min(dayarr)-2454710d
1042                print,'dayarr a',max(dayarr+stash.dayshi)-2454710d,min(dayarr+stash.dayshi)-2454710d
1043      print,'stash.dayshi ',stash.dayshi
1044
1045      print,'fdayarr mx/mn ',max(fdayarr),min(fdayarr)
1046                print,'fdaymin ',fdaymin,' fdaymax ',fdaymax
1047
1048      plot_bad_obs=stash.plot_bad_obs
1049     
1050      typestr=''
1051
1052      if (typeproj eq 2) then begin
1053         xrange=[-999., 999.]
1054         yrange=[0., 90.]
1055      endif
1056      if (typeproj eq 3) then begin
1057         xrange=[-999., 999.]
1058         yrange=[-90, 0.]
1059      endif
1060
1061; add 360 if we're in the pacific
1062
1063      if (xrange(1) gt 180 and typeproj eq 1) then begin
1064                print,'xrange(1) gt 180 ',xrange
1065          wh=where(xarr1 lt 0)
1066                    xarr1(wh)=xarr1(wh)+360
1067                endif
1068
1069      print,'dayarr(0): ',dayarr(0),fdayarr(0), daymin, daymax
1070     
1071;     if (plot_bad_obs eq 1) then begin
1072      wh=where(dep1 ge depmin and dep1 le depmax and $
1073         fdayarr ge fdaymin and fdayarr le fdaymax and $
1074         xarr1 ge xrange(0) and xarr1 le xrange(1) and $
1075         yarr1 ge yrange(0) and yarr1 le yrange(1))
1076;     endif else begin
1077;     wh=where(dep1 ge depmin and dep1 le depmax and $
1078;        fdayarr ge fdaymin and fdayarr le fdaymax and $
1079;        xarr1 ge xrange(0) and xarr1 le xrange(1) and $
1080;        yarr1 ge yrange(0) and yarr1 le yrange(1) and qcarr eq 0)
1081;     endelse     
1082
1083      nelwh=n_elements(wh)
1084
1085      print,'nelwh ',nelwh
1086   
1087      if (wh(0) gt -1) then begin
1088;        obssel=obs1(wh)
1089;        bkgsel=bkg1(wh)
1090                                     
1091         qcsel=qcarr(wh)
1092         lonsel=xarr1(wh)
1093         latsel=yarr1(wh)
1094         daysel=dayarr(wh)
1095         depsel=dep1(wh)
1096         obstypsel=obstypes1(wh)
1097         obnumsel=obnum(wh)
1098         numsel=lonarr(nelwh)
1099
1100         nel=n_elements(obssel)
1101                   
1102; select points here to avoid unnecessary density calculations                   
1103
1104         obssel1=stash.obs(wh)
1105         bkgsel1=stash.bkg(wh)
1106         obssel2=stash.obs2(wh)
1107         bkgsel2=stash.bkg2(wh)
1108
1109         if (density eq 0 and mld eq 0) then begin   
1110            if (salinity eq 0) then begin                 
1111               obssel=obssel1
1112               bkgsel=bkgsel1
1113            endif else begin
1114               obssel=obssel2
1115               bkgsel=bkgsel2
1116            endelse             
1117         endif else begin
1118            if (filetype eq 'CRT') then begin         
1119               bkgsel=sqrt(bkgsel1*bkgsel1 + bkgsel2*bkgsel2)
1120               if (salinity eq 0) then begin
1121                  obssel=sqrt(obssel1*obssel1 + obssel2*obssel2)
1122               endif else begin
1123                  obssel=stash.obs3(wh)               
1124               endelse
1125            endif else begin           ; calculate density if mld flag set
1126               obssel=eos_nemo(obssel1,obssel2)
1127               bkgsel=eos_nemo(bkgsel1,bkgsel2)
1128            endelse
1129         endelse
1130
1131; create an array to store the sum of the square differences for calculating profile rms
1132;        omb2sel=(obssel - bkgsel)^2
1133         nel=n_elements(obssel)
1134      endif else begin        ; if there are no obs then exit subroutine
1135         nel=0
1136         return
1137                endelse 
1138
1139      if ((stash.txt eq 1 or stash.netcdf eq 1) and stash.bindata eq 1) then begin
1140         innovseli=obssel^2
1141                        innovseli2=bkgsel^2
1142                        innovseli3=(obssel-bkgsel)^2
1143                        innovseli4=obssel*bkgsel
1144      endif else begin
1145
1146                             
1147      if (ombtypeplot eq 1) then begin
1148                  typestr='obs - bkg'
1149                        innovseli=obssel-bkgsel
1150                endif
1151      if (ombtypeplot eq 2) then begin
1152                  typestr='obs'
1153                        innovseli=obssel
1154                endif
1155      if (ombtypeplot eq 3) then begin
1156                  typestr='bkg'
1157                        innovseli=bkgsel
1158                endif
1159
1160      if (keyword_set(rmsmean)) then begin
1161                  innovseli=innovseli
1162                        innovseli2=innovseli^2
1163                endif else begin
1164
1165      if (typeplot eq 1) then begin
1166                  typestr='mean '+typestr
1167                endif
1168                if (typeplot eq 2) then begin
1169                  typestr='rms '+typestr
1170                        innovseli=innovseli^2
1171                endif
1172      if (typeplot eq 3) then begin
1173                  typestr='sd '+typestr
1174                        innovseli=innovseli^2      ; ?
1175                endif
1176                if (typeplot eq 4) then begin
1177                  typestr='mean sq '+typestr
1178                        innovseli=innovseli^2
1179                endif
1180
1181      endelse
1182
1183      endelse
1184
1185;      print,'bkgsel ',bkgsel
1186;      print,'obssel ',obssel
1187
1188; tolerance check
1189
1190      if (wh(0) gt -1) then begin
1191         wh2=[-1]
1192
1193;  print,qcsel
1194                   
1195                   if (plot_bad_obs eq 0) then begin
1196;        if (typeplot eq 1 or typeplot eq 2 $
1197;        or typeplot eq 3) then $
1198         if (ombtypeplot eq 1) then $
1199;        wh2=where(bkgsel ne rmdi and obssel ne rmdi)
1200         wh2=where(abs(bkgsel-rmdi) gt abs(rmdi)*0.01 and abs(obssel-rmdi) gt $
1201                           abs(rmdi)*0.01 and qcsel eq 0)   ; tolerance rmdi check
1202         if (ombtypeplot eq 2) then $
1203;                   if (typeplot eq 4) then $
1204;           wh2=where(obssel ne rmdi)
1205         wh2=where(abs(obssel-rmdi) gt abs(rmdi)*0.01 and qcsel eq 0)   ;tol rmdi check
1206         if (ombtypeplot eq 3) then $
1207;        if (typeplot eq 5) then $ 
1208;        wh2=where(bkgsel ne rmdi)
1209         wh2=where(abs(bkgsel-rmdi) gt abs(rmdi)*0.01 and qcsel eq 0)   ;tol rmdi check
1210
1211; check for bad observations
1212;     if (plot_bad_obs eq 1) then begin
1213      endif else if (plot_bad_obs eq 1 and stash.plot_only_bad_obs eq 1) then begin
1214;                 wh2=where(qcsel ne 0 or obssel eq rmdi or bkgsel eq rmdi, count)
1215;                        wh2=where(qcsel ne 0 or abs(obssel-rmdi) lt abs(rmdi)*0.01,count)
1216         wh2=where(qcsel ne 0 and obssel ne rmdi)
1217
1218;                        wh2=where(qcsel ne 0 or abs(obssel-rmdi) lt abs(rmdi)*0.01 or abs(bkgsel-rmdi) lt abs(rmdi)*0.01,count)
1219;                        print,'bad obs: ',count
1220;                        whtmp=where(qcsel eq 0, count)
1221;                        print,'qcsel eq 0: ',count
1222;                        whtmp=where(obssel eq 0, count)
1223;                        print,'obssel eq 0: ',count
1224;                        whtmp=where(bkgsel eq 0, count)
1225;                        print,'bkgsel eq 0: ',count
1226                endif else begin          ; plot all observations
1227;                 wh2=where(obssel)       
1228                        wh2=where(obssel ne rmdi)
1229                endelse
1230
1231         if (wh2(0) gt -1) then begin
1232            obssel=obssel(wh2)
1233            bkgsel=bkgsel(wh2)
1234                      innovseli=innovseli(wh2)
1235;                      if ((stash.txt eq 1 or stash.netcdf eq 1) and stash.bindata eq 1) then begin
1236                      if (n_elements(innovseli2) gt 0) then innovseli2=innovseli2(wh2)
1237                      if (n_elements(innovseli3) gt 0) then innovseli3=innovseli3(wh2)
1238                      if (n_elements(innovseli4) gt 0) then innovseli4=innovseli4(wh2)
1239;           endif
1240            qcsel=qcsel(wh2)
1241            lonsel=lonsel(wh2)
1242            latsel=latsel(wh2)
1243            daysel=daysel(wh2)
1244            depsel=depsel(wh2)
1245            obstypsel=obstypsel(wh2)
1246                      obnumsel=obnumsel(wh2)
1247;                      numsel=numsel(wh2)
1248         endif
1249
1250      print,'nelwh2 ',n_elements(wh2), typeplot, ombtypeplot
1251             
1252      if (wh2(0) gt -1) then begin
1253               
1254                ; calculate vertical gradients if required
1255   if (stash.filetype eq "Prof" and stash.vertgrad eq 1) then $
1256   calcvertgrad, lonsel, latsel, daysel, obssel, bkgsel, qcsel, depsel, rmdi
1257
1258; average profiles on the same latitude and longitude
1259; goes through each point and puts the average in the first instance and flags
1260; the others for later removal
1261; this routine can also be used to bin data
1262      depmx=max(depsel,min=depmn)
1263      print,'depmn/mx: ',depmn, depmx
1264      if (depmx gt depmn or stash.duplicates gt 0 or stash.differences or stash.bindata) then begin
1265
1266; sort by lat and lon
1267; presumably will average in time also...
1268
1269      if (stash.bindata) then begin    ; if we are binning data (in 1 deg bins)
1270                              ; use +0.5 to get nearest round number
1271        latbin=double(fix((latsel+stash.binsize(1)/2.)*(1/stash.binsize(1))))/(1/stash.binsize(1))
1272        lonbin=double(fix((lonsel+stash.binsize(0)/2.)*(1/stash.binsize(0))))/(1/stash.binsize(0))       
1273        latlon=long(latbin*10000000d)+double(lonbin+180d)/360d
1274                endif else begin
1275        latlon=long(double(latsel)*10000000d)+double(lonsel+180d)/360d
1276                endelse
1277               
1278      elsort=sort(latlon)
1279
1280      latsel=latsel(elsort)
1281      lonsel=lonsel(elsort)
1282      obssel=obssel(elsort)
1283                bkgsel=bkgsel(elsort)
1284                innovseli=innovseli(elsort)
1285;                if ((stash.txt eq 1 or stash.netcdf eq 1) and stash.bindata eq 1) then begin
1286                if (n_elements(innovseli2) gt 0) then innovseli2=innovseli2(elsort)
1287                if (n_elements(innovseli3) gt 0) then innovseli3=innovseli3(elsort)
1288                if (n_elements(innovseli4) gt 0) then innovseli4=innovseli4(elsort)
1289;     endif
1290      qcsel=qcsel(elsort)
1291                daysel=daysel(elsort)
1292                depsel=depsel(elsort)
1293      obstypsel=obstypsel(elsort)
1294                obnumsel=obnumsel(elsort)
1295;                numsel=numsel(elsort)
1296      latlon=latlon(elsort)
1297     
1298      nel=n_elements(obssel)
1299                accum=0
1300                isave=0
1301                lon1=-99999.0
1302                lat1=-99999.0
1303                latlon1=-99999.0
1304                day1=-99999
1305                num=0
1306 
1307;     if (nel gt 100) then begin
1308;                info,latsel
1309;     print,latsel(0:100)
1310;                info,lonsel
1311;                print,lonsel(0:100)
1312;                info,latlon
1313;     print,latlon(0:100)
1314 
1315;     latsel(0:100)=latsel(50)
1316;                lonsel(0:100)=lonsel(50)
1317;     latlon=long(double(latsel)*10000000d)+double(lonsel+180d)/360d
1318 
1319;                endif
1320               
1321
1322      print,'DJL ',max(daysel),min(daysel)
1323
1324; loop through the data averaging profiles at the same location               
1325 
1326                sign=1                                           
1327                if (stash.duplicates eq 2) then sign=-1     ; plot the difference rather than the average   
1328 
1329                numsel=lonarr(nel)
1330      for i=0L,nel-1 do begin
1331               
1332                if i mod 100000 eq 0 then print,i,nel
1333
1334;     if (plot_bad_obs eq 0 and obssel(i) ne rmdi) or $
1335;                       (plot_bad_obs eq 1 and lonsel(i) ne rmdi) then begin
1336     
1337;     if (lonsel(i) ne lon1 or latsel(i) ne lat1) then begin
1338
1339; check for new profile or reaching the end of the data
1340; check for the end of a day added
1341
1342      if ((latlon(i) ne latlon1) or (i eq nel-1) or (long(daysel(i)) ne long(day1))) then begin   
1343
1344;     print,i, latlon(i), latlon1, daysel(i), day1               
1345;                if i gt 0 then stop
1346     
1347; average last profile
1348        if (isave ge 0 and num gt 0) then begin
1349               
1350                    if (num gt 1) then begin   
1351               
1352;                 print,num
1353                        numsave=num
1354                        if (stash.duplicates eq 2) then num=1
1355
1356         if (stash.bindata) then begin
1357         latsel(isave)=latsel(isave)/num
1358                        lonsel(isave)=lonsel(isave)/num
1359                        endif                         
1360                  obssel(isave)=obssel(isave)/num
1361                        bkgsel(isave)=bkgsel(isave)/num
1362                        daysel(isave)=daysel(isave)/num
1363                        depsel(isave)=depsel(isave)/num
1364                        innovseli(isave)=innovseli(isave)/num
1365;                        if ((stash.txt eq 1 or stash.netcdf eq 1) and stash.bindata eq 1) then begin
1366                        if (n_elements(innovseli2) gt 0) then innovseli2(isave)=innovseli2(isave)/num
1367                        if (n_elements(innovseli3) gt 0) then innovseli3(isave)=innovseli3(isave)/num
1368                        if (n_elements(innovseli4) gt 0) then innovseli4(isave)=innovseli4(isave)/num
1369;             endif
1370
1371                        num=numsave
1372
1373; save the num
1374         numsel(isave)=num
1375
1376;                        if (fix(lonsel(isave)) eq -7297 and $
1377;                            fix(latsel(isave)) eq -5935) then begin
1378;                         print,'completed ',isave, num                       
1379;         print,'obssel ',obssel(isave)
1380;                         print,'bkgsel ',bkgsel(isave)
1381;                         endif
1382                         
1383                    endif else begin
1384                     numsel(isave)=num
1385                    endelse
1386                    if ((stash.duplicates gt 0 and num eq 1) or (stash.differences and num gt 1)) then begin
1387                        obssel(isave)=rmdi
1388                        bkgsel(isave)=rmdi
1389                        lonsel(isave)=rmdi
1390                        latsel(isave)=rmdi
1391                        qcsel(isave)=rmdi
1392                        innovseli(isave)=rmdi
1393;                        if ((stash.txt eq 1 or stash.netcdf eq 1) and stash.bindata eq 1) then begin
1394                        if (n_elements(innovseli2) gt 0) then innovseli2(isave)=rmdi
1395                        if (n_elements(innovseli3) gt 0) then innovseli3(isave)=rmdi
1396                        if (n_elements(innovseli4) gt 0) then innovseli4(isave)=rmdi
1397;                endif
1398                       
1399                    endif
1400                                           
1401        endif
1402
1403      if (plot_bad_obs eq 0 and obssel(i) ne rmdi) or $
1404                        (plot_bad_obs eq 1 and lonsel(i) ne rmdi) then begin
1405
1406; start a new profile
1407                  lon1=lonsel(i)
1408            lat1=latsel(i)
1409                      latlon1=latlon(i)
1410                      day1=daysel(i)
1411                      isave=i
1412                      accum=1
1413                      num=1
1414
1415      endif    ; plot_bad_obs ...
1416                     
1417 ; accumulating                     
1418                     
1419      endif else begin
1420                 
1421      if (plot_bad_obs eq 0 and obssel(i) ne rmdi) or $
1422                        (plot_bad_obs eq 1 and lonsel(i) ne rmdi) then begin
1423                       
1424                  num=num+1
1425
1426         if (stash.bindata) then begin
1427                        latsel(isave)=latsel(isave)+double(fix((1./stash.binsize(1))*(latsel(i)+stash.binsize(1)/2.)))/(1./stash.binsize(1)) ;latsel(isave)+fix(latsel(i)+0.5)
1428                        lonsel(isave)=lonsel(isave)+double(fix((1./stash.binsize(0))*(lonsel(i)+stash.binsize(0)/2.)))/(1./stash.binsize(0)) ;lonsel(isave)+fix(lonsel(i)+0.5)
1429                        endif
1430
1431                        obssel(isave)=obssel(isave)+obssel(i)*sign
1432                        bkgsel(isave)=bkgsel(isave)+bkgsel(i)*sign
1433                        innovseli(isave)=innovseli(isave)+innovseli(i)*sign
1434;                        if ((stash.txt eq 1 or stash.netcdf eq 1) and stash.bindata eq 1) then begin
1435                        if (n_elements(innovseli2) gt 0) then innovseli2(isave)=innovseli2(isave)+innovseli2(i)*sign
1436                        if (n_elements(innovseli3) gt 0) then innovseli3(isave)=innovseli3(isave)+innovseli3(i)*sign
1437                        if (n_elements(innovseli4) gt 0) then innovseli4(isave)=innovseli4(isave)+innovseli4(i)*sign
1438;             endif
1439
1440                        qcsel(isave)=min([qcsel(isave),qcsel(i)])
1441                        daysel(isave)=daysel(isave)+daysel(i)
1442                        depsel(isave)=depsel(isave)+depsel(i)
1443
1444         if (stash.bindata) then begin
1445                        latsel(i)=rmdi
1446                        lonsel(i)=rmdi               
1447                        endif         
1448                        obssel(i)=rmdi
1449                        bkgsel(i)=rmdi
1450                        innovseli(i)=rmdi
1451;                        if ((stash.txt eq 1 or stash.netcdf eq 1) and stash.bindata eq 1) then begin
1452                        if (n_elements(innovseli2) gt 0) then innovseli2(i)=rmdi
1453                        if (n_elements(innovseli3) gt 0) then innovseli3(i)=rmdi
1454                        if (n_elements(innovseli4) gt 0) then innovseli4(i)=rmdi
1455;             endif
1456                        lonsel(i)=rmdi
1457                        latsel(i)=rmdi
1458                        qcsel(i)=rmdi
1459
1460;                        if (fix(lonsel(isave)*100) eq -7297 and $
1461;                              fix(latsel(isave)*100) eq -5935) then begin
1462;                         print,'isave ',isave, num                       
1463;          print,'obssel ',obssel(isave)
1464;                         print,'bkgsel ',bkgsel(isave)
1465;                        endif
1466
1467      endif    ; plot_bad_obs ...
1468                       
1469                endelse
1470                     
1471;                endif
1472
1473                endfor     ; i  nobs
1474
1475      print,'DJL2x ',max(daysel),min(daysel)
1476
1477
1478      wh5=where(lonsel ne rmdi)           ; clears out but the averaged profiles
1479                if (wh5(0) gt -1) then begin
1480            obssel=obssel(wh5)
1481            bkgsel=bkgsel(wh5)
1482                      innovseli=innovseli(wh5)
1483;                      if ((stash.txt eq 1 or stash.netcdf eq 1) and stash.bindata eq 1) then begin
1484                      if (n_elements(innovseli2) gt 0) then innovseli2=innovseli2(wh5)
1485                      if (n_elements(innovseli3) gt 0) then innovseli3=innovseli3(wh5)
1486                      if (n_elements(innovseli4) gt 0) then innovseli4=innovseli4(wh5)
1487;           endif
1488            qcsel=qcsel(wh5)
1489            lonsel=lonsel(wh5)
1490            latsel=latsel(wh5)
1491            daysel=daysel(wh5)
1492            obstypsel=obstypsel(wh5)
1493                      obnumsel=obnumsel(wh5)
1494                      numsel=numsel(wh5)
1495                      nel=n_elements(obssel)
1496                endif else begin       ; if there are no obs then exit subroutine
1497                  nel=0
1498         return
1499                endelse               
1500      endif
1501
1502      print,'nel wh5 ',n_elements(wh5)
1503      print,'DJL3 ',max(daysel),min(daysel)
1504
1505      typeselect=stash.obstypeselect
1506;                print, 'typeselect ',typeselect
1507      if (typeselect(0) ne "") then begin
1508;                typeselecta=strsplit(typeselect,' ,',/extract)      ; split string on commas or spaces
1509      typeselecta=strsplit(typeselect,',',/extract)      ; split string on commas only
1510                nel=n_elements(typeselecta)
1511                print,'n_elements(typeselecta) ',nel
1512                for i=0L,nel-1 do begin
1513                ; string matching (understands * ?)
1514                st1=strtrim(string(obstypsel),2)      ; strip tailing and leading spaces
1515                st2=strtrim(typeselecta(i),2)
1516                wh6t=where(strmatch(st1,st2,/FOLD_CASE) EQ 1) 
1517;                print, 'wh6t ',typeselecta(i), wh6t
1518      if i eq 0 then wh6=wh6t else wh6=[wh6,wh6t]
1519                endfor
1520;                print,n_elements(wh6)
1521;                print,wh6
1522                wh=where(wh6 gt -1)
1523;                print,wh
1524                if (wh(0) gt -1) then begin
1525                      wh6=wh6(wh)
1526            obssel=obssel(wh6)
1527            bkgsel=bkgsel(wh6)
1528                      innovseli=innovseli(wh6)
1529;                      if ((stash.txt eq 1 or stash.netcdf eq 1) and stash.bindata eq 1) then begin
1530                      if (n_elements(innovseli2) gt 0) then innovseli2=innovseli2(wh6)
1531                      if (n_elements(innovseli3) gt 0) then innovseli3=innovseli3(wh6)
1532                      if (n_elements(innovseli4) gt 0) then innovseli4=innovseli4(wh6)
1533;           endif
1534            qcsel=qcsel(wh6)
1535            lonsel=lonsel(wh6)
1536            latsel=latsel(wh6)
1537            daysel=daysel(wh6)
1538            obstypsel=obstypsel(wh6)
1539                      obnumsel=obnumsel(wh6)
1540                      numsel=numsel(wh6)
1541                      nel=n_elements(obssel)
1542;                      print, 'nel (in obstypeselect) ',nel
1543                endif else begin       ; if there are no obs then exit subroutine
1544                  print,'no matching obs'
1545         nel=0
1546                endelse               
1547      endif                   
1548               
1549                if (stash.printobstypes eq 1) then begin
1550      s=sort(obstypsel)
1551                sobstypsel=obstypsel(s)
1552      u=uniq(sobstypsel)
1553      print,"unique obs types: ",sobstypsel(u)
1554      endif
1555
1556      print,'nel ',nel
1557      stash.symsize=1.0
1558      if (nel gt 250) then stash.symsize=0.75
1559      if (nel gt 1000) then stash.symsize=0.5
1560      if (nel gt 2000) then stash.symsize=0.25
1561      if (nel lt 20) then stash.symsize=2.0
1562
1563      if (!d.name eq 'Z') then stash.symsize=stash.symsize/2
1564
1565      print,'stash.symsize: ',stash.symsize
1566
1567; save data for mean and rms plots
1568; for ombtypeplot=3 typeplot=2 innovsel = (obs - bkg)^2     mean value available from obssel(i) - bkgsel(i)
1569
1570; copy innovseli to innovsel (if we've got some valid data to plot / save)
1571      if (nel gt 0) then innovsel=innovseli
1572
1573      if (stash.txt eq 1 and stash.bindata eq 1) then begin
1574 
1575; The file contains:  longitude, latitude, num of points, obs mean, bkg mean
1576;         mean sq obs,    mean sq bkg,        mean sq obs - bkg,  mean obs * bkg
1577               
1578                outfile=stash.outfile+'_selpoints.txt'
1579                print, 'saving data to ',outfile
1580                OPENW, unit, outfile, /get_lun
1581;                printf,unit,'Output selpoints data: '
1582                printf,unit, 3               ; file version
1583                printf,unit, typeplot
1584                printf,unit, stash.bindata
1585                printf,unit, xrange, yrange
1586                printf,unit, min(daysel), max(daysel)
1587
1588                for i=0L,nel-1 do begin
1589                     printf,unit, lonsel(i), latsel(i), numsel(i), obssel(i), bkgsel(i), $
1590                              innovseli(i), innovseli2(i), innovseli3(i), innovseli4(i), $
1591                              format='(2d10.3, i8, 6d18.8)'
1592                endfor
1593                FREE_LUN, unit
1594               
1595      endif
1596
1597      if (stash.netcdf eq 1 and stash.bindata eq 1) then begin
1598
1599                outfile=stash.outfile+'_selpoints.nc'
1600                print, 'saving data to ',outfile
1601
1602                nid=NCDF_CREATE( outfile, /CLOBBER )
1603               
1604                NCDF_ATTPUT, nid, 'version', 3, /SHORT, /GLOBAL
1605                NCDF_ATTPUT, nid, 'bindata', stash.bindata, /FLOAT, /GLOBAL
1606                NCDF_ATTPUT, nid, 'xrange0', xrange(0), /FLOAT, /GLOBAL
1607                NCDF_ATTPUT, nid, 'xrange1', xrange(1), /FLOAT, /GLOBAL
1608                NCDF_ATTPUT, nid, 'yrange0', yrange(0), /FLOAT, /GLOBAL
1609                NCDF_ATTPUT, nid, 'yrange1', yrange(1), /FLOAT, /GLOBAL
1610                NCDF_ATTPUT, nid, 'mindaysel', min(daysel), /FLOAT, /GLOBAL
1611                NCDF_ATTPUT, nid, 'maxdaysel', max(daysel), /FLOAT, /GLOBAL
1612               
1613                print, 'nel ',nel
1614                nelid = NCDF_DIMDEF( nid, 'n', nel )
1615               
1616                lonid = NCDF_VARDEF ( nid, 'longitudes', [nelid], /FLOAT )
1617                latid = NCDF_VARDEF ( nid, 'latitudes', [nelid], /FLOAT )
1618                numid = NCDF_VARDEF ( nid, 'num', [nelid], /LONG )
1619                obsid = NCDF_VARDEF ( nid, 'obs_mean', [nelid], /FLOAT )
1620                bkgid = NCDF_VARDEF ( nid, 'bkg_mean', [nelid], /FLOAT )
1621                innovseli1id = NCDF_VARDEF ( nid, 'obs_mean_sq', [nelid], /FLOAT )
1622                innovseli2id = NCDF_VARDEF ( nid, 'bkg_mean_sq', [nelid], /FLOAT )
1623                innovseli3id = NCDF_VARDEF ( nid, 'omb_mean_sq', [nelid], /FLOAT )
1624                innovseli4id = NCDF_VARDEF ( nid, 'otimesb_mean', [nelid], /FLOAT )
1625                NCDF_CONTROL, nid, /ENDEF
1626               
1627                if (nel gt 0) then begin
1628                NCDF_VARPUT, nid, lonid, lonsel
1629                NCDF_VARPUT, nid, latid, latsel
1630                NCDF_VARPUT, nid, numid, numsel
1631                NCDF_VARPUT, nid, obsid, obssel
1632                NCDF_VARPUT, nid, bkgid, bkgsel
1633                NCDF_VARPUT, nid, innovseli1id, innovseli
1634                NCDF_VARPUT, nid, innovseli2id, innovseli2
1635                NCDF_VARPUT, nid, innovseli3id, innovseli3
1636                NCDF_VARPUT, nid, innovseli4id, innovseli4
1637                endif
1638               
1639                NCDF_CLOSE, nid
1640               
1641      endif
1642
1643      if (stash.filterout ne '') then begin
1644
1645         ; write out filtered feedback file in new feedback format
1646         ; with some filthy fudges to emulate NEMOVAR format
1647         if (nel gt 0) then begin
1648
1649            print, 'writing filtered feedback file ',stash.filterout
1650
1651            ; split type string in two to be used as STATION_IDENTIFIER and STATION_TYPE
1652            typesel=strarr(nel)
1653            instsel=strarr(nel)
1654
1655            for iob = 0, nel-1 do begin
1656               tempsel=strsplit(obstypsel[iob],/extract)
1657               typesel[iob]=tempsel[0]
1658               ; catch for data without space and instrument id in types (e.g. altimeter)
1659               instsel[iob]=''
1660               if (n_elements(tempsel) gt 1 ) then $
1661                  instsel[iob]=tempsel[1]
1662            endfor
1663
1664            ; setup netcdf file
1665            nid=NCDF_CREATE( stash.filterout, /CLOBBER )
1666
1667            ; global attribute to emulate NEMOVAR feedback files
1668            NCDF_ATTPUT, nid, 'title', "NEMO observation operator output", /GLOBAL
1669            if (stash.obstypeselect ne '' ) then $
1670               NCDF_ATTPUT, nid, 'filter', stash.obstypeselect, /GLOBAL
1671            NCDF_ATTPUT, nid, 'minday', min(daysel), /FLOAT, /GLOBAL
1672            NCDF_ATTPUT, nid, 'maxday', max(daysel), /FLOAT, /GLOBAL
1673
1674            nobid = NCDF_DIMDEF( nid, 'N_OBS', nel )
1675            ndepid = NCDF_DIMDEF( nid, 'N_LEVELS', 1 ) ; temporarily input 1d profile arrays
1676
1677            if (stash.filetype eq 'feedback') then begin
1678
1679               ncvars = stash.varname
1680               nnvars = n_elements(ncvars)
1681
1682            endif else begin
1683
1684               nnvars = 1
1685               ncvars = stash.filetype
1686               if (stash.filetype eq 'Prof') then begin
1687                  ncvars = 'POTM'
1688                  if (salinity eq 1) then ncvars = 'PSAL'
1689                  nextid = NCDF_DIMDEF( nid, 'N_EXTRA', 1 )
1690               endif
1691
1692            endelse
1693
1694            nvarid = NCDF_DIMDEF( nid, 'N_VARS', nnvars )
1695
1696            nentid = NCDF_DIMDEF( nid, 'N_ENTRIES', 1 )
1697
1698            nstrnamid = NCDF_DIMDEF( nid, 'STRINGNAM', 8 )
1699            nstrwmoid = NCDF_DIMDEF( nid, 'STRINGWMO', 8 )
1700            nstrtypid = NCDF_DIMDEF( nid, 'STRINGTYP', 4 )
1701            nstrjulid = NCDF_DIMDEF( nid, 'STRINGJULD', 14 )
1702
1703            varid = NCDF_VARDEF ( nid, 'VARIABLES', [nstrnamid,nvarid], /CHAR )
1704            entid = NCDF_VARDEF ( nid, 'ENTRIES', [nstrnamid,nentid], /CHAR )
1705            signid = NCDF_VARDEF ( nid, 'STATION_IDENTIFIER', [nstrwmoid,nobid], /CHAR )
1706            typeid = NCDF_VARDEF ( nid, 'STATION_TYPE', [nstrtypid,nobid], /CHAR )
1707
1708            lonid = NCDF_VARDEF ( nid, 'LONGITUDE', [nobid], /DOUBLE )
1709            latid = NCDF_VARDEF ( nid, 'LATITUDE', [nobid], /DOUBLE )
1710            depid = NCDF_VARDEF ( nid, 'DEPTH', [nobid], /DOUBLE )
1711;            depid = NCDF_VARDEF ( nid, 'DEPTH', [ndepid,nobid], /DOUBLE )
1712            julid = NCDF_VARDEF ( nid, 'JULD', [nobid], /DOUBLE )
1713            jrefid = NCDF_VARDEF ( nid, 'JULD_REFERENCE', [nstrjulid], /CHAR )
1714
1715            obsqcid = NCDF_VARDEF ( nid, 'OBSERVATION_QC', [nobid], /LONG )
1716
1717            FillVal=9999
1718            NCDF_ATTPUT, nid, depid, '_Fillvalue', FillVal, /DOUBLE
1719            NCDF_ATTPUT, nid, obsqcid, '_Fillvalue', FillVal, /LONG
1720
1721            obsids=intarr(nnvars)
1722            bkgids=intarr(nnvars)
1723            qcids=intarr(nnvars)
1724            lqcids=intarr(nnvars)
1725
1726            for ivar = 0, nnvars-1 do begin
1727
1728               obsids[ivar] = NCDF_VARDEF ( nid, ncvars[ivar]+'_OBS', [nobid], /FLOAT )
1729               NCDF_ATTPUT, nid, obsids[ivar], '_Fillvalue', FillVal, /FLOAT
1730               bkgids[ivar] = NCDF_VARDEF ( nid, ncvars[ivar]+'_Hx', [nobid], /FLOAT )
1731               NCDF_ATTPUT, nid, bkgids[ivar], '_Fillvalue', FillVal, /FLOAT
1732               qcids[ivar] = NCDF_VARDEF ( nid, ncvars[ivar]+'_QC', [nobid], /LONG )
1733               NCDF_ATTPUT, nid, qcids[ivar], '_Fillvalue', FillVal, /FLOAT
1734               lqcids[ivar] = NCDF_VARDEF ( nid, ncvars[ivar]+'_LEVEL_QC', [nobid], /LONG )
1735               NCDF_ATTPUT, nid, lqcids[ivar], '_Fillvalue', FillVal, /FLOAT
1736
1737            endfor
1738
1739            NCDF_CONTROL, nid, /ENDEF
1740
1741            JulDRef = '19500101000000'
1742            NCDF_VARPUT, nid, jrefid, JulDRef
1743            RefDate = JULDAY(fix(strmid(JulDRef,4,2)), fix(strmid(JulDRef,6,2)), fix(strmid(JulDRef,0,4)), $
1744                             fix(strmid(JulDRef,8,2)), fix(strmid(JulDRef,10,2)), fix(strmid(JulDRef,12,2)))
1745            NCDF_VARPUT, nid, julid,daysel-RefDate
1746
1747            NCDF_VARPUT, nid, varid, ncvars
1748            NCDF_VARPUT, nid, entid, 'Hx'
1749
1750            NCDF_VARPUT, nid, lonid, lonsel
1751            NCDF_VARPUT, nid, latid, latsel
1752
1753            NCDF_VARPUT, nid, depid, fltarr(nel)
1754            NCDF_VARPUT, nid, obsqcid, fltarr(nel)+1
1755
1756            for ivar = 0, nnvars-1 do begin
1757
1758               NCDF_VARPUT, nid, obsids[ivar], obssel
1759               NCDF_VARPUT, nid, bkgids[ivar], bkgsel
1760               NCDF_VARPUT, nid, qcids[ivar], qcsel+1
1761               NCDF_VARPUT, nid, lqcids[ivar], fltarr(nel)+1
1762
1763            endfor
1764
1765            NCDF_VARPUT, nid, typeid, typesel
1766            NCDF_VARPUT, nid, signid, instsel
1767
1768         endif
1769
1770         NCDF_CLOSE, nid
1771
1772      endif ; filterout
1773
1774      if (nel gt 0) then begin   
1775
1776         if (typeplot eq 2) then begin
1777                        innovsel=sqrt(innovsel)
1778                   endif
1779         if (typeplot eq 3) then begin
1780;  std dev = sqrt ( mean x2 - (mean x)^2)
1781         x2=innovsel
1782         print,'typeplot ',typeplot,' ombtypeplot ',ombtypeplot
1783         if (ombtypeplot eq 1) then $
1784                           innovsel=sqrt(x2-(obssel-bkgsel)^2)
1785                        if (ombtypeplot eq 2) then $
1786                           innovsel=sqrt(x2-(obssel)^2)                       
1787                        if (ombtypeplot eq 3) then $
1788                           innovsel=sqrt(x2-(bkgsel)^2)                       
1789                   endif
1790
1791         if (keyword_set(rmsmean)) then begin
1792                     innovsel=innovseli                  ; mean
1793                        if (n_elements(innovseli2) gt 0) then innovsel2=innovseli2  ; rms                 
1794                   endif
1795
1796
1797         innovmean=total(innovsel)/nel
1798         innovmean2=total(innovsel^2)/nel
1799         innovsd=sqrt(innovmean2-innovmean^2)
1800         print,'mean innovsel ',innovmean
1801         print,'sd innovsel ',innovsd
1802
1803      endif
1804
1805     
1806
1807      endif
1808
1809      endif
1810
1811      wh=where(numsel eq 0)
1812                if (wh(0) gt -1) then numsel(wh)=1
1813
1814END
1815
1816PRO setupct, r, g, b, coltable=coltable, white_on_black=white_on_black
1817
1818if (coltable eq 0 or coltable lt -2) then begin
1819; get color table and modify
1820                loadct, 13
1821                stretch, -40, 256
1822;                tvlct,r,g,b,/get               
1823;                r(0)=0
1824;                g(0)=0
1825;                b(0)=0
1826;                r(255)=255
1827;                g(255)=255
1828;                b(255)=255             
1829;                tvlct,r,g,b
1830endif else begin
1831   if (coltable eq -1) then begin
1832      restore_colors, 'spectrum.clr'
1833   endif else if (coltable eq -2) then begin
1834                restore_colors, '~frbe/spectrum_alt.xdr',/asis
1835        endif else begin
1836      loadct,coltable
1837        endelse           
1838endelse
1839
1840tvlct,r,g,b,/get
1841i0=0
1842i1=255
1843if (n_elements(white_on_black) eq 1 and !d.name ne "PS") then begin
1844   if (white_on_black eq 0) then begin
1845         i0=255
1846                i1=0   
1847        endif
1848endif
1849r(i0)=0
1850g(i0)=0
1851b(i0)=0
1852r(i1)=255
1853g(i1)=255
1854b(i1)=255             
1855tvlct,r,g,b
1856
1857END
1858
1859PRO plotpoints, stash
1860
1861      nplots=1
1862                if (stash.pmulti eq 2) then nplots=stash.pmulti
1863
1864      for nplot=0,nplots-1 do begin
1865
1866      if (stash.busy eq 1) then return
1867      stash.busy = 1
1868
1869      num_cols=254
1870
1871      print,"!d.name ",!d.name
1872      IF (stash.drawID gt -1) then begin
1873      IF (!d.name eq 'X' or !d.name eq 'Z') then begin
1874                        device, decomposed=0
1875           endif
1876                IF (!d.name eq 'X') then WSET, stash.drawID
1877      endif
1878
1879      IF (!d.name ne 'Z') then setupct, r, g, b, $
1880                   coltable=stash.coltable,white_on_black=stash.white_on_black      ; setup color table
1881
1882      noerase=0
1883      if (stash.pmulti eq 2) then begin
1884                  if (stash.pmultinum gt 0) then noerase=1
1885                  !p.multi=[stash.pmultinum,2,1]
1886                        stash.pmultinum=(stash.pmultinum+1) mod 2
1887                endif else begin
1888                  !p.multi=0
1889                endelse
1890               
1891;     xarr1=stash.xarr
1892;     yarr1=stash.yarr
1893;     dep1=stash.dep
1894;     obs1=stash.obs
1895;     bkg1=stash.bkg
1896;     depmin=stash.depmin
1897      xrange=stash.xrange
1898      yrange=stash.yrange
1899;     dayarr=stash.dayarr
1900      daymin=stash.daymin
1901      daymax=stash.daymax
1902;     typeplot=stash.typeplot
1903      typeproj=stash.typeproj
1904;                print,"stash.map_file: ",stash.map_file
1905;     if (strlen(stash.map_file) gt 0) then map_file=stash.map_file
1906                               
1907;     print,depmin,xrange,yrange
1908
1909      dummy=[0,0]
1910     
1911     
1912;     plot,dummy,/nodata,yrange=[-90,90],xrange=[-180,180]
1913;     plot,dummy,/nodata,yrange=[-50,80],xrange=[-110,40],$
1914;        xstyle=1,ystyle=1, xtitle='Longitude',ytitle='Latitude'
1915
1916      if (typeproj ne 1) then begin
1917            origlon=0.
1918         if (typeproj eq 2) then begin
1919                           origlat=90.
1920                                if (yrange(0) lt 0) then yrange(0)=45
1921                                if (yrange(1) lt 0) then yrange(1)=90
1922                        endif
1923         if (typeproj eq 3) then begin
1924                           origlat=-90.
1925                                if (yrange(0) gt 0) then yrange(0)=-90
1926                                if (yrange(1) gt 0) then yrange(1)=-45                               
1927                        endif     
1928         ;scale=0.35
1929                        smap = map_proj_init('Polar Stereographic')
1930      endif
1931
1932      ; select points to plot
1933      typestr=""
1934               
1935spawn,'echo plotpoints 1 `date`'
1936               
1937      selpoints, stash, lonsel, latsel, innovsel, qcsel, daysel, obstypsel, obnumsel, numsel, typestr
1938               
1939spawn,'echo plotpoints 2 `date`'
1940
1941      if (stash.txt eq 0 and stash.netcdf eq 0) then begin
1942      nelsel=n_elements(innovsel)
1943
1944                nptsstr=strtrim(string(nelsel),2)
1945                print,nptsstr, nelsel
1946
1947;     if (nelsel gt 0) then begin
1948
1949
1950
1951      ;convert lats and lons to projection positions
1952;                if (nelsel gt 0) then begin
1953;                print,'min/max lonsel: ',min(lonsel),max(lonsel)
1954;     if (typeproj ne 1) then begin
1955;                  coords = map_proj_forward(lonsel,latsel,MAP=smap)
1956;       lonsel=coords(0,*)
1957;       latsel=coords(1,*)
1958;     endif
1959;                endif
1960
1961
1962
1963      print,'nelsel: ',nelsel
1964      nelinnov=n_elements(innovsel)
1965      print,'nelinnov: ',nelinnov
1966
1967                jul_to_dtstr, daymin, dateminstr, /notime
1968                jul_to_dtstr, daymax, datemaxstr, /notime     
1969
1970                strtsal=''
1971;                print,'stash.filetype ',stash.filetype
1972                if (stash.filetype eq "Prof" or stash.filetype eq "feedback") then begin
1973                     strtsal='T:  '
1974                     if (stash.salinity eq 1) then strtsal='S:  '
1975                     if (stash.density  eq 1) then strtsal='Density:  '
1976                     if (stash.mld eq 1) then strtsal='MLD:  '
1977                endif else begin
1978                     if (stash.filetype eq "CRT") then begin
1979                        strtsal='U:  '
1980                        if (stash.salinity eq 1) then strtsal='V:  '
1981                        if (stash.density  eq 1) then strtsal='Speed:  '
1982                        if (stash.salinity eq 1 and stash.density  eq 1) then strtsal='Total Speed:  '
1983                     endif else begin
1984                        strtsal=stash.filetype+':  '
1985                     endelse
1986                endelse
1987                if (stash.vertgrad eq 1) then strtsal='Grad '+strtsal
1988      nptsstr2='Points: '+nptsstr+'  '
1989      titlestr=strtsal+typestr+': '+dateminstr+' to '+datemaxstr
1990      print,'titlestr ',titlestr
1991
1992; get max and min values
1993      mx=0
1994                mn=0
1995                meaninnov=0
1996                rmsinnov=0
1997      if (nelinnov gt 0) then begin
1998; NB these values are just the rms and mean of the points plotted...
1999; and do not take account of the datapoints used in profiles plotted...
2000         mx=max(innovsel)
2001                   mn=min(innovsel)
2002         meaninnov=avg(innovsel)
2003                   rmsinnov=sqrt(avg(innovsel^2))
2004;                   n=total(numsel)
2005;                   wh=where(numsel gt 0)
2006;                   print, 'innovsel: ',max(innovsel), min(innovsel)
2007;                   print, 'numsel: ',max(numsel), min(numsel)
2008;                   if (n gt 0) then begin
2009;                      x=total(innovsel(wh)*numsel(wh))
2010;                      x2=total(innovsel(wh)^2*numsel(wh))
2011;                      meaninnov=x/n
2012;                      rmsinnov=sqrt(x2/n)
2013;                   endif
2014                endif
2015
2016      subtitle=''
2017      subtitle=nptsstr2+'depths: '+strtrim(string(long(stash.depmin)),2)+'-'+$
2018                  strtrim(string(long(stash.depmax)),2)
2019                if (stash.obstypeselect ne "") then subtitle=subtitle+'  filtered type: '+stash.obstypeselect       
2020      subtitle=subtitle+'  extrema: '+$
2021                  strtrim(string(mn,format='(G0.4)'),2)+', '+strtrim(string(mx,format='(G0.4)'),2)
2022      subtitle=subtitle+'  mean: '+strtrim(string(meaninnov,format='(G0.4)'))+$
2023                  ' rms: '+strtrim(string(rmsinnov,format='(G0.4)'))
2024     
2025;                print,'noerase ',noerase
2026               
2027      if (typeproj ne 1) then begin                 
2028                 map_set, /STEREOGRAPHIC, origlat, 0, /continents, title=titlestr+'!C'+subtitle, $
2029                         ymargin=[10.,5.], /label, latlab=xrange(0), lonlab=yrange(0), $
2030                         limit=[yrange(0),xrange(0),yrange(1),xrange(1)],/isotropic, hires=stash.hires_map, noerase=noerase
2031                endif else begin       
2032                  P0lon=0
2033                  P0lat=0
2034                  if (xrange(1) gt 180) then P0lon=180
2035                  print,'ranges: ',yrange(0),xrange(0),yrange(1),xrange(1)
2036                  map_set, P0lat, P0lon, /continents, title=titlestr+'!C'+subtitle, $
2037                         ymargin=[10.,5.], /label, latlab=xrange(0), lonlab=yrange(0), $
2038                         limit=[yrange(0),xrange(0),yrange(1),xrange(1)],/isotropic, hires=stash.hires_map, noerase=noerase 
2039                endelse
2040
2041                                                                 
2042      if (nelsel eq 0 or nelinnov eq 0) then begin
2043;        pp_contour,fld,title=titlestr,/nodata, /proportional, map_file=map_file
2044      endif else begin
2045;        pp_contour,fld,title=titlestr,/nodata, /proportional, map_file=map_file
2046                         
2047                if (stash.fmx ne stash.rmdi) then mx=stash.fmx
2048                if (stash.fmn ne stash.rmdi) then mn=stash.fmn
2049
2050      print,'setting stash mx/mn'
2051      stash.mx=mx
2052                stash.mn=mn
2053
2054;     levs=contour_levels([mn,mx],nlevels=15)
2055;     levs=findgen(10+1)/10.*(mx-mn)+mn
2056
2057      print, 'mn mx ',mn,mx
2058
2059
2060
2061;; De-select points outside window?
2062;     
2063;    xnd= [ -!X.s(0), 1. ]/!X.s(1)
2064;    ynd= [ -!Y.s(0), 1. ]/!Y.s(1)
2065;    ov=0                ; overlap
2066;    clip_data= [ (!ppp.position(0)-ov)*xnd(1) + xnd(0) , $
2067;                 (!ppp.position(1)-ov)*ynd(1) + ynd(0) , $
2068;                 (!ppp.position(2)+ov)*xnd(1) + xnd(0) , $
2069;                 (!ppp.position(3)+ov)*ynd(1) + ynd(0) ]
2070;
2071;  wh=where(lonsel ge clip_data(0) and lonsel le clip_data(2)  and $
2072;               latsel ge clip_data(1) and latsel  le clip_data(3))
2073;         if (wh(0) gt -1) then begin
2074;           latsel=latsel(wh)
2075;     lonsel=lonsel(wh)
2076;     innovsel=innovsel(wh)
2077;         endif
2078
2079
2080      print,'daymin=',daymin
2081      print,'daymax=',daymax
2082;     print,dayarr(wh)
2083
2084   
2085;     titlestr=typestr+': '+titlestr
2086     
2087
2088
2089      mx=max(innovsel)
2090                mn=min(innovsel)
2091                if (stash.fmx ne stash.rmdi) then mx=stash.fmx
2092                if (stash.fmn ne stash.rmdi) then mn=stash.fmn
2093               
2094;                stash.mx=mx
2095;                stash.mn=mn
2096               
2097;     levs=contour_levels([mn,mx],nlevels=30)
2098      levs=findgen(15)/15.*(mx-mn)+mn
2099                               
2100                print, 'typeproj: ',typeproj
2101     
2102;     wh=sort(innovsel)
2103;     print,'innovsel ',innovsel(wh)
2104
2105      clr=fix((innovsel-mn)/(mx-mn)*(num_cols-1))+1
2106
2107; prevent colours going out of range
2108      wh=where(clr lt 1)
2109                if (wh(0) gt -1) then clr(wh)=1
2110                wh=where(clr gt num_cols)
2111                if (wh(0) gt -1) then clr(wh)=num_cols
2112
2113;     print,'clr ',clr(wh)
2114
2115                case stash.sym of
2116      1: USERSYM, [-1,-1,1,1,-1],[1,-1,-1,1,1], /FILL
2117                2: usersym, cos(2.0*!pi*indgen(17)/16.0), sin(2.0*!pi*indgen(17)/16.0), thick=2.0, /fill
2118                endcase
2119                 
2120                 spawn,'echo plotpoints 3 `date`'
2121                 
2122;;    rgb=transpose(reform(color24_to_rgb(clr),nelwh,3))
2123
2124      symsize=stash.symsize*stash.symscale 
2125;     for i=0L,nelsel-1 do $
2126;          plots,lonsel(i),latsel(i),psym=8,color=clr,symsize=symsize
2127
2128      if (stash.bindata eq 0) then begin
2129         plots,lonsel,latsel,psym=8,color=clr,symsize=symsize,noclip=0
2130;                plots,lonsel,latsel,psym=3,color=clr
2131
2132;     print, 'lonsel latsel'
2133;     print,lonsel, format='(4g20.10)'
2134;                print,latsel, format='(4g20.10)'
2135      endif
2136               
2137               
2138;                print, 'lonseld latseld'
2139;     print,lonsel - lonsel(0), format='(4g20.10)'
2140;                print,latsel - latsel(0), format='(4g20.10)'
2141
2142;     wh=where(lonsel eq lonsel(0) and latsel eq latsel(0))
2143;                print,wh
2144
2145      endelse
2146
2147                 spawn,'echo plotpoints 4 `date`'
2148      print,"stash.symsize ",stash.symsize
2149
2150      print,'limit=',[yrange(0),xrange(0),yrange(1),xrange(1)]
2151
2152;     endif       ; nel lonsel > 0
2153               
2154      if (stash.bindata) then begin
2155                     fillval=-32000
2156                     mxlat=max(latsel)
2157                     mnlat=min(latsel)
2158                     mxlon=max(lonsel)
2159                     mnlon=min(lonsel)
2160                     dlon=stash.binsize(0)
2161                     dlat=stash.binsize(1)
2162                     nlats=(mxlat-mnlat)/dlat+1
2163                     nlons=(mxlon-mnlon)/dlon+1
2164                     print,'nlons ',nlons,' nlats ',nlats
2165                     if (nlats ge 2 and nlons ge 2) then begin
2166                     innovsel2=fltarr(nlons,nlats)
2167;                     latsel2=fltarr(nlons,nlats)
2168;                     lonsel2=fltarr(nlons,nlats)
2169         latsel2=findgen(nlats)*dlat+mnlat
2170                        lonsel2=findgen(nlons)*dlon+mnlon
2171                     innovsel2(*)=fillval       ; missing data
2172;                     latsel2(*)=fillval
2173;                     lonsel2(*)=fillval
2174                     nelin=n_elements(innovsel)
2175                     print,mnlat,mnlon
2176                     iiarr=(lonsel-mnlon)/dlon
2177                     ijarr=(latsel-mnlat)/dlat
2178                     print,min(ijarr),max(ijarr),min(ijarr),max(ijarr)
2179                     info,latsel2
2180                     info,lonsel2
2181                     for i=0L,nelin-1 do begin
2182                     ii=iiarr(i)
2183                     iip=ii+1
2184                     if (iip ge nlons-1) then iip=ii
2185                     iim=ii-1
2186                     if (iim lt 0) then iim=0
2187                     ij=ijarr(i)
2188                     ijp=ij+1
2189                     if (ijp ge nlats-1) then ijp=ij
2190                     ijm=ij-1
2191                     if (ijm lt 0) then ijm=0
2192
2193                     innovsel2(ii,ij)=innovsel(i)
2194                     ; fill in gaps
2195                     if (innovsel2(iip,ij) eq fillval) then innovsel2(iip,ij)=innovsel(i)
2196                     if (innovsel2(ii,ijp) eq fillval) then innovsel2(ii,ijp)=innovsel(i)
2197                     if (innovsel2(iim,ij) eq fillval) then innovsel2(iim,ij)=innovsel(i)
2198                     if (innovsel2(ii,ijm) eq fillval) then innovsel2(ii,ijm)=innovsel(i)
2199;                     latsel2(ii,ij)=latsel(i)
2200;                     lonsel2(ii,ij)=lonsel(i)
2201                     endfor
2202                     info,latsel2
2203                     print,latsel2
2204                     info,lonsel2
2205                     print,lonsel2
2206                     
2207                     TVLCT, R, G, B, /get
2208           TVLCT,R(1),G(1),B(1),0        ; bodge fix for wrong black contour   
2209
2210         wh=where(innovsel2 gt mx-(mx-mn)*0.001 and innovsel2 ne fillval)
2211                        if (wh(0) gt -1) then innovsel2(wh)=mx-(mx-mn)*0.001
2212         wh=where(innovsel2 lt mn and innovsel2 ne fillval)
2213                        if (wh(0) gt -1) then innovsel2(wh)=mn
2214
2215           levs2=nice_contourlevels([mn,mx],nlevels=50)
2216
2217           contour,innovsel2,lonsel2,latsel2,levels=levs2, /overplot, /cell_fill, min_value=-1000
2218                     
2219                     TVLCT,R(0),G(0),B(0),0     ; restore original colour table
2220                     endif
2221                     
2222      endif
2223
2224      if (daymin eq daymax and daymin eq 0) then begin
2225                  print,'No data file'
2226                        xyouts, 0.15, 0.5, 'No data files', /normal, charsize=8
2227                endif
2228
2229; plot colorbar
2230      levs=nice_contourlevels([mn,mx],nlevels=10)
2231
2232      print, levs
2233               
2234      barclr=findgen(num_cols)+1
2235                colorbar_idl, [mn,mx], barclr, levs
2236      endif          ; stash.txt eq 0
2237
2238      stash.busy=0
2239
2240                if (stash.pmulti eq 2) then begin
2241                  typeproj=stash.typeproj
2242                        print, 'old stash.typeproj: ',stash.typeproj
2243                  if (typeproj eq 2) then stash.typeproj=3     ;plot north and south poles
2244                  if (typeproj eq 3) then stash.typeproj=2
2245                        print, 'new stash.typeproj: ',stash.typeproj
2246                endif
2247
2248      endfor         ; nplots
2249
2250
2251end
2252
2253;------------------------------------------
2254;Profile plotting window and event handling
2255;------------------------------------------
2256
2257PRO profilewindow_event, ev
2258  WIDGET_CONTROL, ev.TOP, GET_UVALUE=stash2
2259  WIDGET_CONTROL, ev.ID, GET_UVALUE=uval
2260   print,uval
2261   if (uval eq "PRINTPW") then begin
2262      ps=1
2263      eps=0
2264      landscape=1
2265      pr2,file="~/plotprofile.ps",landscape=landscape,ps=ps,eps=eps,color=1
2266      plotprofile,stash2
2267      prend2,/view
2268   endif
2269        if (uval eq "SAVEPW") then begin
2270      thisDevice = !D.Name
2271         psave=!p
2272
2273         Set_Plot, 'Z'        ; do graphics in the background
2274      Device, Set_Resolution=[640,480]    ; clear any existing stuff
2275         !p.charsize=0.75
2276
2277      if (stash2.white_on_black eq 0) then begin
2278; flip background and foreground color
2279                        pcolor=!p.color
2280                        pbackground=!p.background
2281                        !p.color=pbackground
2282                        !p.background=pcolor
2283
2284         endif
2285         print,'!p.color,!p.background ',!p.color,!p.background
2286
2287         setupct, r, g, b, $
2288                   coltable=stash2.coltable,white_on_black=stash2.white_on_black    ; setup color table
2289; plot data
2290      plotprofile,stash2
2291      snapshot = TVRD()
2292         WRITE_GIF,"~/plotprofile.gif",snapshot, r, g, b
2293         Device, Z_Buffer=1      ; reset graphics mode
2294         Set_Plot, thisDevice
2295      !p=psave
2296       
2297         spawn,'xv ~/plotprofile.gif'
2298       
2299        endif
2300        if (uval eq "TXTPW") then begin
2301         plotprofile,stash2,/txt
2302        endif
2303end
2304
2305PRO plotprofile,stash2, txt=txt
2306
2307   if (n_elements(txt) eq 0) then txt=0
2308
2309;** need to deal with multiple times
2310
2311   dep2=stash2.dep2
2312   val2=stash2.val2
2313   obs2=stash2.obs2
2314   bkg2=stash2.bkg2
2315   obstype2=stash2.obstype2
2316        qcarr2=stash2.qcarr2
2317        dayarr2=stash2.dayarr2
2318   datestr=stash2.datestr
2319   xselstr=stash2.xselstr
2320   yselstr=stash2.yselstr
2321   rmdi=stash2.rmdi
2322
2323        multstr=""
2324   mindayarr2=min(dayarr2,max=maxdayarr2)
2325       
2326        whs=sort(dayarr2)
2327        u=uniq(dayarr2)
2328
2329        nu=n_elements(u)
2330        dayarr2u=dayarr2(u)
2331
2332        print,'dayarr2u: ',dayarr2u
2333       
2334        if (maxdayarr2 ne mindayarr2) then multstr="mult times "+strtrim(string(nu),2)
2335
2336;  print,'obstype2 ',obstype2
2337
2338   print,'profile plot_bad_obs: ',stash2.plot_bad_obs
2339   print, obs2
2340        print, bkg2
2341        print, qcarr2
2342
2343      if (stash2.plot_bad_obs eq 0) then begin
2344      wh=where(obs2 ne rmdi and qcarr2 eq 0)
2345                endif else begin
2346;                wh=where(obs2 eq obs2)
2347      wh=where(obs2 ne rmdi)
2348                endelse
2349               
2350      if (wh(0) gt -1) then begin
2351      depo1=dep2(wh)
2352      valo1=val2(wh)
2353      obso1=obs2(wh)
2354      bkgo1=bkg2(wh)
2355                qcarro1=qcarr2(wh)
2356      endif
2357               
2358               
2359      if (stash2.plot_bad_obs eq 0) then begin
2360      wh=where(bkg2 ne rmdi and qcarr2 eq 0)
2361      endif else begin
2362;                wh=where(bkg2 eq bkg2)
2363      wh=where(bkg2 ne rmdi)
2364                endelse
2365
2366      if (wh(0) gt -1) then begin
2367      depb1=dep2(wh)
2368      valb1=val2(wh)
2369      obsb1=obs2(wh)
2370      bkgb1=bkg2(wh)
2371                qcarrb1=qcarr2(wh)
2372      endif
2373               
2374
2375;loop thru times
2376
2377      if (n_elements(bkgb1) gt 0) then begin
2378                  xmx=max([obso1,bkgb1],min=xmn)
2379                        ymx=max([depo1,depb1],min=ymn)
2380                endif else begin
2381                  xmx=max(obso1,min=xmn)
2382                        ymx=max(depo1,min=ymn)
2383                endelse
2384
2385; set max and min depth based on selection of depths from the main window
2386      if (ymx gt stash2.depmax) then ymx=stash2.depmax
2387                if (ymn lt stash2.depmin) then ymn=stash2.depmin
2388
2389; reset xrange based on new ymx and ymn
2390      if (n_elements(bkgb1) gt 0) then begin
2391                        wh1=where(depo1 ge ymn and depo1 le ymx)
2392                        wh2=where(depb1 ge ymn and depb1 le ymx)
2393                  if (wh1(0) gt -1 and wh2(0) gt -1) then xmx=max([obso1(wh1),bkgb1(wh2)],min=xmn)
2394                endif else begin
2395                        wh1=where(depo1 ge ymn and depo1 le ymx)
2396                  if (wh1(0) gt -1) then xmx=max(obso1(wh1),min=xmn)
2397                endelse
2398
2399; add a little bit to xrange to make plots look better
2400      if (xmn eq xmx) then begin
2401                  xmn=xmn-0.5
2402                        xmx=xmx+0.5
2403                endif
2404      dxrange=(xmx-xmn)*0.02
2405                xmx=xmx+dxrange
2406                xmn=xmn-dxrange
2407
2408; add a little bit to yrange to make plots look better
2409      if (ymn eq ymx) then begin
2410                  ymn=ymn-0.5
2411                        ymx=ymx+0.5
2412                endif
2413      dyrange=(ymx-ymn)*0.01
2414                ymx=ymx+dyrange
2415                ymn=ymn-dyrange
2416      print,'** yrange: ',ymn, ymx
2417
2418      if (txt eq 1) then begin
2419              outfile='dataplot_profile.txt'
2420         print, 'saving data to ',outfile
2421         OPENW, unit, outfile, /get_lun
2422                endif
2423   
2424      for t=0,nu-1 do begin               ; loop through times
2425
2426      wht=where(dayarr2 eq dayarr2u(t))
2427      depo=depo1(wht)
2428                valo=valo1(wht)
2429                obso=obso1(wht)
2430                bkgo=bkgo1(wht)
2431                qcarro=qcarro1(wht)
2432                dayarro=dayarr2(wht)
2433      if (n_elements(depb1) gt 0) then depb=depb1(wht)
2434                if (n_elements(valb1) gt 0) then valb=valb1(wht)
2435                if (n_elements(obsb1) gt 0) then obsb=obsb1(wht)
2436                if (n_elements(bkgb1) gt 0) then bkgb=bkgb1(wht)
2437                if (n_elements(qcarrb1) gt 0) then qcarrb=qcarrb1(wht)
2438               
2439
2440
2441;print data
2442
2443;     if (txt eq 1) then begin
2444;     print,'rmdi: ',rmdi
2445;                print,'profile data'
2446;     for i=0,n_elements(depo)-1 do begin
2447;                 print,depo(i),obso(i),bkgo(i),qcarro(i)
2448;                endfor
2449;                print,'profile data'
2450;     for i=0,n_elements(depb)-1 do begin
2451;                 print,depb(i),obsb(i),bkgb(i),qcarrb(i)
2452;                endfor
2453;     endif
2454
2455      typestr="type: "+string(obstype2(0))
2456   
2457;get variable type
2458      vartype=size(obstype2,/type)
2459      if (vartype ne 7) then begin
2460         if (obstype2(0) eq rmdi) then typestr=""
2461      endif
2462
2463; line thickness
2464      thick1=2
2465                thick2=2
2466
2467;     color1='FFCC66'x
2468      color1=102
2469      color2=!p.color
2470                color3=234
2471               
2472                linestyle1=0
2473                linestyle2=0
2474               
2475      who=sort(depo)
2476                whb=sort(depb)
2477
2478      obso_srto=obso(who)
2479                depo_srto=depo(who)
2480                qcarro_srto=qcarro(who)
2481                whbad=where(qcarro ne 0)
2482               
2483                if (txt eq 0) then begin
2484               
2485                if (t eq 0) then $
2486      plot,obso_srto,depo_srto,ystyle=1,xstyle=1,linestyle=linestyle1,$
2487         xrange=[xmn,xmx],yrange=[ymx,ymn],thick=thick1, $
2488         psym=-4,ytitle='Level',xtitle='Value', $
2489         title=datestr+" "+multstr+" ("+xselstr+","+yselstr+")   "+typestr, /nodata
2490               
2491                oplot,obso_srto,depo_srto,linestyle=linestyle1,psym=-4, thick=thick1, color=color1
2492                if (whbad(0) gt -1) then begin
2493                  oplot,obso_srto(whbad),depo_srto(whbad),psym=4, thick=thick1, color=color3
2494                        oplot,obso_srto(whbad),depo_srto(whbad),psym=1, thick=thick1, color=color3
2495                endif
2496      if (n_elements(bkgb) gt 0) then oplot,bkgb(whb),depb(whb),psym=-5,thick=thick2, color=color2,$
2497                  linestyle=linestyle2
2498      xcoord=0.8
2499      ycoord=0.2
2500;                if (stash2.salinity eq 1) then xcoord=0.15 ; left corner legend for salinity
2501      nel=n_elements(who)
2502                if (obso(who(0)) lt obso(who(nel-1))) then xcoord=0.15
2503                if (t eq 0) then $
2504
2505                  ycoord2=ycoord-0.05
2506                  xcoord2=xcoord+0.03
2507                  xcoord3=xcoord+0.05
2508                  plots, [xcoord,xcoord2],[ycoord,ycoord], psym=-4, linestyle=linestyle1, /normal, $
2509                     thick=thick1, color=color1
2510                  xyouts, xcoord+0.05, ycoord, 'obs', /normal
2511                  plots, [xcoord,xcoord2],[ycoord2,ycoord2], psym=-5, linestyle=linestyle2, /normal, $
2512                     thick=thick2, color=color2
2513                  xyouts, xcoord+0.05, ycoord2, 'bkg',/normal               
2514               
2515                endif else begin ; txt
2516                 
2517                  printf, unit,'observations ',n_elements(who)             
2518                  for i=0L,n_elements(who)-1 do begin             
2519                     printf, unit, depo(i), obso(i), dayarro(i), format='(3f18.5)'
2520                  endfor
2521                  printf, unit, 'background ',n_elements(who)
2522                  for i=0L,n_elements(whb)-1 do begin
2523                     printf, unit, depb(i), bkgb(i), dayarro(i), format='(3f18.5)'
2524                  endfor             
2525                               
2526                endelse
2527                                       
2528                endfor  ; t
2529               
2530                if (txt eq 1) then FREE_LUN, unit
2531               
2532end
2533
2534PRO profilewindow,dep2,val2,obs2,bkg2,obstype2,qcarr2,dayarr2,datestr,xselstr,yselstr,rmdi,salinity=salinity, $
2535      plot_bad_obs=plot_bad_obs, density=density, white_on_black=white_on_black, coltable=coltable, $
2536                depmax=depmax, depmin=depmin
2537   basepw=WIDGET_BASE(/column)
2538   drawpw = WIDGET_DRAW(basepw, XSIZE=640, YSIZE=480)
2539        buttonBase = Widget_Base(basepw, Row=1)
2540   buttonpw = WIDGET_BUTTON(buttonBase, VALUE='Print',uvalue="PRINTPW")
2541   button2pw = WIDGET_BUTTON(buttonBase, VALUE='Save',uvalue="SAVEPW")
2542        button3pw = WIDGET_BUTTON(buttonBase, VALUE='Text file',uvalue="TXTPW")
2543   WIDGET_CONTROL, basepw, /REALIZE
2544   
2545        if (n_elements(salinity) eq 0) then salinity=0
2546        if (n_elements(density) eq 0) then density=0
2547   if (n_elements(plot_bad_obs) eq 0) then plot_bad_obs=0
2548        if (n_elements(white_on_black) eq 0) then white_on_black=1
2549
2550      stash2 = { dep2:dep2,val2:val2,obs2:obs2, bkg2:bkg2, $
2551         obstype2:obstype2, qcarr2:qcarr2, dayarr2:dayarr2, $
2552         datestr:datestr,xselstr:xselstr, yselstr:yselstr, $
2553         rmdi:rmdi, salinity:salinity, plot_bad_obs:plot_bad_obs,$
2554                        density:density, white_on_black:white_on_black, $
2555                        coltable:coltable, depmax:depmax, depmin:depmin }
2556
2557      plotprofile,stash2
2558
2559   WIDGET_CONTROL, basepw, SET_UVALUE=stash2
2560
2561   XMANAGER, 'profilewindow', basepw, /NO_BLOCK
2562end
2563
2564;---------------------------------------
2565;Worst points, window and event handling
2566;---------------------------------------
2567
2568PRO worstpoints,stash2
2569
2570;stash2 is a combination of main stash and worstpoints specific stuff
2571
2572        xrange=stash2.xrange
2573        yrange=stash2.yrange
2574        selpoints,stash2,lonsel, latsel, innovsel, qcsel, daysel, obstypsel, obnumsel, obnumsel, numsel, typestr
2575
2576   str1='Worst points in '+$
2577        '   lons ('+strtrim(string(xrange(0)),2)+','+strtrim(string(xrange(1)),2)+$
2578                  ')   lats ('+strtrim(string(yrange(0)),2)+','+strtrim(string(yrange(1)),2)+')'+$
2579        '   depths '+strtrim(string(long(stash2.depmin)),2)+'-'+strtrim(string(long(stash2.depmax)),2)               
2580        print,str1
2581
2582   WIDGET_CONTROL, stash2.wwlabel1, set_value=str1
2583
2584        count=n_elements(innovsel)
2585        wh=sort(abs(innovsel))
2586        wh0=reverse(wh)       ; reverse order starting with the largest       
2587        cmax=min([count,10])
2588   stash2.lonselw(0:cmax-1)=lonsel(wh0(0:cmax-1))
2589        stash2.latselw(0:cmax-1)=latsel(wh0(0:cmax-1))
2590        stash2.dayselw(0:cmax-1)=daysel(wh0(0:cmax-1))
2591        print,'lon      lat       '+typestr+'     qc         date    '
2592        for j1=1,cmax do begin
2593           j=j1-1
2594       jul_to_dtstr,daysel(wh0(j)),datestr
2595;           datedt=jul_to_dt(daysel(wh0(j)))
2596;           datestr=strtrim(fix(datedt.year),2)+'/'+ $
2597;       strtrim(string(fix(datedt.month),format='(i2.2)'),2)+'/'+ $
2598;       strtrim(string(fix(datedt.day),format='(i2.2)'),2)+' '+ $
2599;       strtrim(string(fix(datedt.hour),format='(i2.2)'),2)+':'+ $
2600;       strtrim(string(fix(datedt.minute),format='(i2.2)'),2)+':'+ $
2601;       strtrim(string(fix(datedt.second),format='(i2.2)'),2)
2602     
2603           print,lonsel(wh0(j)),latsel(wh0(j)),innovsel(wh0(j)),$
2604                qcsel(wh0(j)),' ',datestr
2605
2606        endfor
2607       
2608        if (n_elements(stash2) gt 0) then begin
2609        i=0
2610   WIDGET_CONTROL, stash2.labels2(i,0), set_value="Lon" & i=i+1
2611   WIDGET_CONTROL, stash2.labels2(i,0), set_value="Lat" & i=i+1
2612   WIDGET_CONTROL, stash2.labels2(i,0), set_value=typestr & i=i+1
2613   WIDGET_CONTROL, stash2.labels2(i,0), set_value="QC" & i=i+1
2614        WIDGET_CONTROL, stash2.labels2(i,0), set_value="Date" & i=i+1
2615       
2616        for j1=1,cmax do begin
2617           j=j1-1
2618           jul_to_dtstr,daysel(wh0(j)),datestr
2619;           datedt=jul_to_dt(daysel(wh0(j)))
2620;           datestr=strtrim(fix(datedt.year),2)+'/'+ $
2621;       strtrim(string(fix(datedt.month),format='(i2.2)'),2)+'/'+ $
2622;       strtrim(string(fix(datedt.day),format='(i2.2)'),2)+' '+ $
2623;       strtrim(string(fix(datedt.hour),format='(i2.2)'),2)+':'+ $
2624;       strtrim(string(fix(datedt.minute),format='(i2.2)'),2)+':'+ $
2625;       strtrim(string(fix(datedt.second),format='(i2.2)'),2)
2626
2627        i=0
2628   WIDGET_CONTROL, stash2.labels2(i,j1), $
2629         set_value=strtrim(string(lonsel(wh0(j)),format='(f9.3)'),2) & i=i+1
2630   WIDGET_CONTROL, stash2.labels2(i,j1), $
2631         set_value=strtrim(string(latsel(wh0(j)),format='(f9.3)'),2) & i=i+1
2632   WIDGET_CONTROL, stash2.labels2(i,j1), set_value=strtrim(string(innovsel(wh0(j))),2) & i=i+1
2633   WIDGET_CONTROL, stash2.labels2(i,j1), set_value=strtrim(string(qcsel(wh0(j))),2) & i=i+1
2634        WIDGET_CONTROL, stash2.labels2(i,j1), set_value=strtrim(datestr,2) & i=i+1
2635
2636;           print,lonsel(wh0(j)),latsel(wh0(j)),innovsel(wh0(j)),$
2637;                qcsel(wh0(j)),' ',datestr
2638
2639        endfor
2640
2641       
2642        endif
2643           
2644
2645
2646
2647end
2648
2649pro worstpointswindow,stash
2650   basepw=WIDGET_BASE(/column)
2651;  drawpw = WIDGET_DRAW(basepw, XSIZE=640, YSIZE=480)
2652
2653   wwlabel1 = WIDGET_LABEL(basepw, XSIZE=480, VALUE="Worst points")
2654               
2655        nj=11
2656        ni=5
2657   labels=intarr(nj)
2658        for j=0,nj-1 do begin       
2659          labels(j)=Widget_Base(basepw,row=1)
2660        endfor
2661
2662        labsiz=[50,50,50,50,125]
2663   labels2=intarr(ni+2,nj)
2664        for j=0,nj-1 do begin
2665         for i=0,ni-1 do begin
2666         labels2(i,j) = WIDGET_LABEL(labels(j), XSIZE=labsiz(i), $
2667    VALUE="a")
2668    endfor
2669         if (j gt 0) then begin
2670         labels2(i,j) = WIDGET_BUTTON(labels(j), VALUE='Zoom to', uvalue="Zoom to "+strtrim(string(j),2))
2671         labels2(i+1,j) = WIDGET_BUTTON(labels(j), VALUE='Profile', uvalue="Profile "+strtrim(string(j),2))
2672         endif
2673        endfor
2674
2675        stashb={ wwlabel1:wwlabel1, labels2:labels2, ni:ni, nj:nj, lonselw:dblarr(10), latselw:dblarr(10), dayselw:dblarr(10)}
2676   stash2=CREATE_STRUCT(stash,stashb)
2677
2678   buttonpw = WIDGET_BUTTON(basepw, VALUE='Print',uvalue="PRINTPW")
2679               
2680   WIDGET_CONTROL, basepw, /REALIZE
2681        worstpoints,stash2
2682   WIDGET_CONTROL, basepw, SET_UVALUE=stash2
2683   XMANAGER, 'worstpointswindow', basepw, /NO_BLOCK
2684end
2685
2686PRO worstpointswindow_event, ev
2687  WIDGET_CONTROL, ev.TOP, GET_UVALUE=stash2
2688  WIDGET_CONTROL, ev.ID, GET_UVALUE=uval
2689   print,uval
2690
2691      xarr1=stash2.xarr
2692      yarr1=stash2.yarr
2693      dep1=stash2.dep
2694      dayarr=stash2.dayarr
2695      daymin=stash2.daymin
2696      daymax=stash2.daymax
2697                obstypes1=stash2.obstypes
2698                xrange=stash2.xrange
2699      yrange=stash2.yrange
2700      xs=xrange(1)-xrange(0)
2701      ys=yrange(1)-yrange(0)
2702      rmdi=stash2.rmdi
2703                if (stash2.salinity eq 0) then begin
2704        obs1=stash2.obs
2705        bkg1=stash2.bkg
2706                  qcarr1=stash2.qcarr
2707                endif else begin
2708        obs1=stash2.obs2
2709        bkg1=stash2.bkg2
2710                  qcarr1=stash2.qcarr2
2711      endelse                 
2712
2713   if (uval eq "PRINTTSW") then begin
2714      ps=1
2715      eps=0
2716      landscape=1
2717      pr2,file="~/worstpoints.ps",landscape=landscape,ps=ps,eps=eps,color=1
2718      worstpoints,stash2
2719      prend2,/view
2720   endif
2721        print,strmid(uval,0,7)
2722        if (strcmp(strmid(uval,0,7),"Zoom to")) then begin
2723         num=fix(strmid(uval,7,3))
2724                xsel=stash2.lonselw(num-1)
2725                ysel=stash2.latselw(num-1)
2726                print,num,xsel,ysel 
2727                dataplot_zoom,stash2, xsel, ysel, 1
2728                plotpoints,stash2             
2729                           
2730        endif
2731        if (strcmp(strmid(uval,0,7),"Profile")) then begin
2732         num=fix(strmid(uval,7,3))
2733                xsel=stash2.lonselw(num-1)
2734                ysel=stash2.latselw(num-1)
2735                daysel=stash2.dayselw(num-1)
2736;                print,num,xsel,ysel,daysel,daysel-long(daysel),float(daysel)
2737;                info,dayarr
2738;                info,daysel
2739;     info,xarr1
2740;                info,ysel
2741;                info,yarr1
2742;                info,xsel
2743               
2744                wh1=where(xarr1 eq xsel and yarr1 eq ysel)
2745;                print,'dayarr(wh1) ',dayarr(wh1), float(dayarr(wh1))
2746;                print,'dayarr(wh1)-long ',dayarr(wh1)-long(dayarr(wh1))
2747;                print,'dep1(wh1) ',dep1(wh1)
2748                wh2=where(xarr1 eq xsel and yarr1 eq ysel and float(dayarr) eq float(daysel))
2749;                print,'wh2 ',wh2
2750;                if (wh2(0) gt -1) then print,'daysel(wh2) ',daysel(wh2)
2751       
2752         wh=where(xarr1 eq xsel and yarr1 eq ysel and long(dayarr) eq long(daysel))
2753                if (wh(0) gt -1) then begin
2754      dep2=dep1(wh)
2755      val2=abs(obs1(wh)-bkg1(wh))
2756      obs2=obs1(wh)
2757      bkg2=bkg1(wh)
2758      obstype2=obstypes1(wh)
2759                qcarr2=qcarr1(wh)
2760                dayarr2=dayarr(wh)
2761      jul_to_dtstr,stash2.dayselw(num-1),datestr
2762      xselstr=string(xsel)
2763      yselstr=string(ysel)
2764
2765      profilewindow,dep2,val2,obs2,bkg2,obstype2,qcarr2,dayarr2,datestr,xselstr,yselstr,rmdi, salinity=stash2.salinity, $
2766                  white_on_black=stash.white_on_black, coltable=stash.coltable, depmax=stash.depmax, depmin=stash.depmin
2767                endif
2768   endif
2769
2770; store values
2771  WIDGET_CONTROL, ev.TOP, SET_UVALUE=stash2
2772
2773end
2774
2775PRO filterwindow, stash
2776
2777   val=stash.obstypeselect
2778
2779   obstypes=stash.obstypes
2780        uniquetypes=strjoin(obstypes(uniq(obstypes, sort(obstypes))),string(10b))
2781       
2782   base = WIDGET_BASE(GROUP_LEADER=stash.base ,/modal,/column)
2783       
2784       
2785   intextid = CW_FIELD(base, TITLE = "Filter type", /FRAME, value=val, uvalue='intext', /RETURN_EVENTS)
2786        outtextid = WIDGET_TEXT(base, value="Choose from:"+string(10b)+uniquetypes, $
2787         scr_xsize=200, scr_ysize=200, /scroll)
2788
2789   buttonBase = Widget_Base(base, Row=1)
2790      cancelID = Widget_Button(buttonBase, Value='Cancel', uvalue='cancel')
2791      acceptID = Widget_Button(buttonBase, Value='Accept', uvalue='accept')
2792
2793   WIDGET_CONTROL, base, /REALIZE
2794
2795   ptr = Ptr_New({text:"", nocancel:0})
2796
2797   stash2={ ptr:ptr, intextid:intextid }
2798        WIDGET_CONTROL, base, SET_UVALUE=stash2, /No_Copy
2799
2800   XMANAGER, 'filterwindow', base       
2801
2802        theText = (*ptr).text
2803        nocancel = (*ptr).nocancel
2804        Ptr_Free, ptr
2805       
2806        print, 'finished ',theText, nocancel
2807
2808; if obstypeselect has changed then replot the points
2809               
2810        if (stash.obstypeselect ne theText and nocancel eq 1) then begin
2811                stash.obstypeselect=theText
2812         plotpoints,stash
2813        endif
2814       
2815       
2816end
2817
2818PRO filterwindow_event, ev
2819  WIDGET_CONTROL, ev.TOP, GET_UVALUE=stash2
2820  WIDGET_CONTROL, ev.ID, GET_UVALUE=uval
2821  WIDGET_CONTROL, ev.ID, GET_VALUE=val
2822;     info,uval
2823;        info,val
2824   print,uval
2825        print,val
2826
2827   if (uval eq "accept" or uval eq "intext") then begin
2828;        stash2.val=val
2829
2830   Widget_Control, stash2.intextid, Get_Value=theText
2831   print,'the text ',theText
2832
2833        (*stash2.ptr).text=theText
2834        (*stash2.ptr).nocancel=1
2835        endif         
2836       
2837;   WIDGET_CONTROL, ev.TOP, SET_UVALUE=stash2       
2838       
2839        WIDGET_CONTROL, ev.top, /DESTROY
2840
2841end
2842
2843;--------------------------------------
2844; Input area
2845;--------------------------------------
2846
2847PRO inputareawindow, xrange, yrange, toplevel, success
2848
2849        val1=xrange(0)
2850        val2=xrange(1)
2851        val3=yrange(0)
2852        val4=yrange(1)
2853       
2854   base = WIDGET_BASE(GROUP_LEADER=toplevel,/modal,/column)
2855               
2856   intextid1 = CW_FIELD(base, TITLE = "lon1", /FRAME, value=val1, uvalue='intext1', /RETURN_EVENTS)
2857   intextid2 = CW_FIELD(base, TITLE = "lon2", /FRAME, value=val2, uvalue='intext2', /RETURN_EVENTS)
2858   intextid3 = CW_FIELD(base, TITLE = "lat1", /FRAME, value=val3, uvalue='intext3', /RETURN_EVENTS)
2859   intextid4 = CW_FIELD(base, TITLE = "lat2", /FRAME, value=val4, uvalue='intext4', /RETURN_EVENTS)
2860
2861;       outtextid = WIDGET_TEXT(base, value="Choose from:"+string(10b)+uniquetypes, $
2862;           scr_xsize=200, scr_ysize=200, /scroll)
2863
2864   buttonBase = Widget_Base(base, Row=1)
2865      acceptID = Widget_Button(buttonBase, Value='Accept', uvalue='accept')
2866      cancelID = Widget_Button(buttonBase, Value='Cancel', uvalue='cancel')
2867
2868   WIDGET_CONTROL, base, /REALIZE
2869
2870   ptr = Ptr_New({text1:"", text2:"", text3:"", text4:"", nocancel:0})
2871
2872   stash2={ ptr:ptr, intextid1:intextid1,  intextid2:intextid2, $
2873            intextid3:intextid3, intextid4:intextid4}
2874        WIDGET_CONTROL, base, SET_UVALUE=stash2, /No_Copy
2875
2876   XMANAGER, 'inputareawindow', base       
2877
2878        lon1 = (*ptr).text1
2879        lon2 = (*ptr).text2
2880        lat1 = (*ptr).text3
2881        lat2 = (*ptr).text4
2882        nocancel = (*ptr).nocancel
2883        Ptr_Free, ptr
2884       
2885        print, 'finished ',lon1, lon2, lat1, lat2, nocancel
2886
2887; if obstypeselect has changed then replot the points
2888               
2889        if (nocancel eq 1) then begin
2890      xrange=[float(lon1),float(lon2)]
2891                yrange=[float(lat1),float(lat2)]
2892                xrange=xrange(sort(xrange))
2893                yrange=yrange(sort(yrange))
2894                wh=where(xrange gt 360)
2895                if (wh(0) gt -1) then xrange(wh)=360
2896                wh=where(xrange lt -180)
2897                if (wh(0) gt -1) then xrange(wh)=-180
2898                wh=where(yrange gt 90)
2899                if (wh(0) gt -1) then yrange(wh)=90
2900                wh=where(yrange lt -90)
2901                if (wh(0) gt -1) then yrange(wh)=-90
2902               
2903                if (xrange(0) eq xrange(1)) then xrange(1)=xrange(1)+0.1
2904                if (yrange(0) eq yrange(1)) then yrange(1)=yrange(1)+0.1
2905               
2906                success=1
2907        endif
2908       
2909       
2910end
2911
2912PRO inputareawindow_event, ev
2913  WIDGET_CONTROL, ev.TOP, GET_UVALUE=stash2
2914  WIDGET_CONTROL, ev.ID, GET_UVALUE=uval
2915  WIDGET_CONTROL, ev.ID, GET_VALUE=val
2916;     info,uval
2917;        info,val
2918   print,uval
2919        print,val
2920
2921   if (uval eq 'accept') then begin
2922;        stash2.val=val
2923
2924   Widget_Control, stash2.intextid1, Get_Value=theText1
2925   Widget_Control, stash2.intextid2, Get_Value=theText2
2926        Widget_Control, stash2.intextid3, Get_Value=theText3
2927        Widget_Control, stash2.intextid4, Get_Value=theText4
2928   print,'the text ',theText1, theText2, theText3, theText4
2929
2930        (*stash2.ptr).text1=theText1
2931        (*stash2.ptr).text2=theText2
2932        (*stash2.ptr).text3=theText3
2933        (*stash2.ptr).text4=theText4
2934        (*stash2.ptr).nocancel=1
2935        endif         
2936       
2937;   WIDGET_CONTROL, ev.TOP, SET_UVALUE=stash2       
2938       
2939        WIDGET_CONTROL, ev.top, /DESTROY
2940
2941end
2942
2943;--------------------------------------
2944; Input max min
2945;--------------------------------------
2946
2947PRO inputmxmnwindow, fmx, fmn, mx, mn, rmdi, toplevel, success
2948
2949   print,'fmx/mn ',fmx,fmn
2950        print,'mx/mn ',mx,mn
2951
2952        val1=mx
2953        val2=mn
2954
2955   base = WIDGET_BASE(GROUP_LEADER=toplevel,/modal,/column)
2956               
2957   intextid1 = CW_FIELD(base, TITLE = "mx", /FRAME, value=val1, uvalue='intext1', /RETURN_EVENTS)
2958   intextid2 = CW_FIELD(base, TITLE = "mn", /FRAME, value=val2, uvalue='intext2', /RETURN_EVENTS)
2959
2960   text="Max and min not locked"
2961        if (fmx ne rmdi) then text="Max locked"
2962        if (fmn ne rmdi) then text="Min Locked"
2963        if (fmx ne rmdi and fmn ne rmdi) then text="Max and min locked"
2964        outtextid = WIDGET_TEXT(base, value=text, $
2965         scr_xsize=200, scr_ysize=40)
2966
2967   buttonBase = Widget_Base(base, Row=1)
2968      acceptID = Widget_Button(buttonBase, Value='Accept/lock', uvalue='accept')
2969           resetID = Widget_Button(buttonBase, Value='Reset/free', uvalue='reset')
2970      cancelID = Widget_Button(buttonBase, Value='Cancel', uvalue='cancel')
2971           
2972
2973   WIDGET_CONTROL, base, /REALIZE
2974
2975   ptr = Ptr_New({text1:"", text2:"", nocancel:0})
2976
2977   stash2={ ptr:ptr, intextid1:intextid1,  intextid2:intextid2}
2978        WIDGET_CONTROL, base, SET_UVALUE=stash2, /No_Copy
2979
2980   XMANAGER, 'inputmxmnwindow', base       
2981
2982        fmx = float((*ptr).text1)
2983        fmn = float((*ptr).text2)
2984
2985   info, fmx
2986        info, fmn
2987
2988        nocancel = (*ptr).nocancel
2989        Ptr_Free, ptr
2990       
2991        print, 'finished ',fmx,fmn, mx, mn
2992
2993; if obstypeselect has changed then replot the points
2994               
2995        if (nocancel eq 1) then begin               
2996                success=1
2997        endif
2998   if (nocancel eq 2) then begin
2999         fmx=rmdi
3000                fmn=rmdi
3001                success=1
3002        endif
3003 
3004end
3005
3006PRO inputmxmnwindow_event, ev
3007  WIDGET_CONTROL, ev.TOP, GET_UVALUE=stash2
3008  WIDGET_CONTROL, ev.ID, GET_UVALUE=uval
3009  WIDGET_CONTROL, ev.ID, GET_VALUE=val
3010;     info,uval
3011;        info,val
3012   print,'uval ',uval,' val ',val
3013
3014   if (uval eq 'accept') then begin
3015;        stash2.val=val
3016
3017   Widget_Control, stash2.intextid1, Get_Value=theText1
3018   Widget_Control, stash2.intextid2, Get_Value=theText2
3019   print,'the text ',theText1, theText2
3020
3021        (*stash2.ptr).text1=theText1
3022        (*stash2.ptr).text2=theText2
3023        (*stash2.ptr).nocancel=1
3024        endif         
3025
3026   if (uval eq 'reset') then begin
3027        (*stash2.ptr).nocancel=2
3028        endif
3029       
3030;   WIDGET_CONTROL, ev.TOP, SET_UVALUE=stash2       
3031       
3032        WIDGET_CONTROL, ev.top, /DESTROY
3033
3034end
3035
3036
3037
3038pro infowindow
3039   basepw=WIDGET_BASE(/column)
3040   xsz=360
3041
3042        iwlabel0 = WIDGET_LABEL(basepw, XSIZE=xsz, VALUE="Dataplot")
3043        iwlabela = WIDGET_LABEL(basepw, XSIZE=xsz, VALUE="--------")
3044        iwlabelb = WIDGET_LABEL(basepw, XSIZE=xsz, VALUE="")     
3045   iwlabel1 = WIDGET_LABEL(basepw, XSIZE=xsz, VALUE="To report bugs or for help")
3046        iwlabel2 = WIDGET_LABEL(basepw, XSIZE=xsz, VALUE="contact Dan Lea")
3047        iwlabel3 = WIDGET_LABEL(basepw, XSIZE=xsz, VALUE="daniel.lea@metoffice.gov.uk")
3048                 
3049   WIDGET_CONTROL, basepw, /REALIZE
3050
3051;  WIDGET_CONTROL, basepw, SET_UVALUE=stash2
3052;  XMANAGER, 'worstpointswindow', basepw, /NO_BLOCK
3053end
3054
3055PRO dataplot_zoom, stash, xdata, ydata, zoominout
3056
3057      xrange=stash.xrange
3058      yrange=stash.yrange
3059      xrangedef=stash.xrangedef
3060      yrangedef=stash.yrangedef
3061   
3062
3063      print,xrange,yrange
3064
3065;     zx=xrange(0)
3066;     zy=yrange(0)
3067      sx=xrange(1)-xrange(0)
3068      sy=yrange(1)-yrange(0)
3069      sxd=xrangedef(1)-xrangedef(0)
3070      syd=yrangedef(1)-yrangedef(0)
3071
3072      print,sx,sy
3073
3074      if (zoominout eq 1) then begin
3075      sx=sx/2.
3076      sy=sy/2.
3077; try to squarify for zoomed in version
3078      sx2=(sx+sy)/2.
3079;     sx=max(sx,sx2)
3080;     sy=max(sy,sx2)
3081;     sx=(sx+sx2)/2.
3082;     sy=(sy+sx2)/2.
3083
3084; try to tend to initial proportions for zooming in
3085      sx=(sx+1.0*(sxd*0.5/(sxd+syd)))/2.
3086                sy=(sy+1.0*(syd*0.5/(sxd+syd)))/2.
3087
3088      endif
3089      if (zoominout eq -1) then begin
3090      sx=sx*2.
3091      sy=sy*2.
3092 ; try to make similar to initial proportions for zooming out               
3093                sx2=(sx+sy)*2.
3094      sx=(sx+sx2*(sxd*0.5/(sxd+syd)))/2.
3095      sy=(sy+sx2*(syd*0.5/(sxd+syd)))/2.
3096               
3097      endif
3098     
3099
3100      print,'sx sy ',sx,sy, xdata, ydata
3101
3102                oxrange=xrange               
3103     
3104      xrange=[xdata-sx/2.,xdata+sx/2.]
3105      yrange=[ydata-sy/2.,ydata+sy/2.]
3106
3107      if (oxrange(1) le 180) then begin      ; 0 longitude centred
3108      if (xrange(0) lt xrangedef(0)) then xrange(0)=xrangedef(0)
3109      if (xrange(1) gt xrangedef(1)) then xrange(1)=xrangedef(1)
3110      endif else begin           ; 180 longitude centred
3111      if (xrange(0) lt xrangedef(0)+180) then xrange(0)=xrangedef(0)+180
3112      if (xrange(1) gt xrangedef(1)+180) then xrange(1)=xrangedef(1)+180
3113      endelse               
3114
3115      if (yrange(0) lt yrangedef(0)) then yrange(0)=yrangedef(0)
3116      if (yrange(1) gt yrangedef(1)) then yrange(1)=yrangedef(1)
3117
3118      print,'xrange ',xrange,' yrange ',yrange
3119               
3120                stash.xrange=xrange
3121      stash.yrange=yrange
3122
3123end
3124
3125;---------------------------------------
3126;Time series plotting and event handling
3127;---------------------------------------
3128
3129PRO plottimeseries, stash
3130 
3131      xrange=stash.xrange
3132      yrange=stash.yrange
3133;     dayarr=stash.dayarr
3134      daymin=stash.daymin
3135      daymax=stash.daymax
3136      dayminl=stash.dayminl
3137      daymaxl=stash.daymaxl
3138
3139      numtimeseries=stash.numtimeseries      ; if 1 then plot number of points
3140
3141;     print,depmin,xrange,yrange
3142
3143;     wh=where(dep1 eq ev.index)
3144;     wh=where(dep1 eq depmin and $
3145;        dayarr ge daymin and dayarr le daymax)
3146
3147
3148
3149;     nelwh=n_elements(wh)
3150
3151      ; select points to plot
3152      typestr=''
3153
3154      if (stash.plottssel eq 1) then begin
3155         dayminl=daymin
3156         daymaxl=daymax
3157      endif
3158
3159      print,'djl ',dayminl,daymaxl, min(stash.dayarr), max(stash.dayarr)
3160
3161      selpoints, stash, lonsel, latsel, innovsel, qcsel, daysel, obstypsel, obnumsel, numsel, typestr, $
3162                  daymin=dayminl, daymax=daymaxl, /rmsmean, innovsel2=innovsel2
3163
3164      print,'djl ',min(daysel), max(daysel)
3165
3166      if (stash.plottssel eq 1) then begin
3167                  daymaxl=daymaxl+1
3168                endif
3169         
3170
3171;     plot,daysel, ystyle=1
3172
3173      print,'n_elements(daysel) ',n_elements(daysel), min(daysel), max(daysel)
3174      print, dayminl, daymaxl
3175
3176; group in 1/4 day bins
3177
3178;     binspday=double(4.0)    ; every 6 hours
3179;                binspday=double(8.0)    ; every 3 hours
3180      binspday=double(stash.binspday)
3181      nbins=(daymaxl-dayminl)*binspday+1
3182
3183      print, daymaxl, dayminl, binspday
3184
3185                x2arr=dblarr(nbins)
3186                xarr=dblarr(nbins)
3187                narr=dblarr(nbins)
3188                tarr=dblarr(nbins)
3189      meants=dblarr(nbins)
3190                rmsts=dblarr(nbins)
3191 
3192      print, 'nbins: ',nbins
3193
3194      if n_elements(innovsel) gt 0 then begin
3195
3196      print,max(daysel),min(daysel)
3197 
3198                for ibin=0L,nbins-1 do begin
3199                  timmn=double(dayminl)+ibin/binspday
3200                  timmx=timmn+1/binspday
3201;                  print, ibin, timmn, timmx, timmx-timmn
3202                  wh=where(daysel ge timmn and daysel lt timmx, count)
3203                  print, 'ibin ',ibin, timmn, timmx, count
3204; redo based on the number of observations
3205                  if (wh(0) ne -1) then begin
3206                  innovsels=innovsel(wh)
3207                  innovsel2s=innovsel2(wh)
3208                  numsels=numsel(wh)
3209                 
3210        x2arr(ibin)=total(innovsel2s*numsels)
3211        xarr(ibin)=total(innovsels*numsels)
3212;       narr(ibin)=n_elements(innovsels)
3213        narr(ibin)=total(numsels)
3214;                  print,ibin, x2arr(ibin), xarr(ibin)
3215                  endif
3216                  tarr(ibin)=timmn
3217                endfor
3218
3219      meants=xarr/narr
3220                rmsts=sqrt(x2arr/narr)
3221
3222;     print,'meants: ',meants
3223;     print,'rmsts: ',rmsts
3224
3225      endif ; n_elements(innovsel)
3226
3227      wh=where(finite(meants))      ; find finite values
3228      ymx=max([meants(wh),rmsts(wh)],min=ymn)
3229
3230;put a bit of white space around ymx,ymn
3231      dymxmn=(ymx-ymn)*0.02
3232                ymx=ymx+dymxmn
3233                ymn=ymn-dymxmn
3234
3235      print,'ymx/mn ',ymx,ymn
3236
3237                strtsal=''
3238;                print,'stash.filetype ',stash.filetype
3239                if (stash.filetype eq 'Prof' or stash.filetype eq 'feedback') then begin
3240                  strtsal='T:  '
3241                        if (stash.salinity eq 1) then strtsal='S:  '
3242                endif else begin
3243                  strtsal=stash.filetype+':  '
3244                endelse
3245
3246                title=strtsal+typestr+'   lons ('+string(xrange(0))+','+string(xrange(1))+$
3247                  ')   lats ('+string(yrange(0))+','+string(yrange(1))+')' 
3248                print,title           
3249
3250;                dtarr=jul_to_dt(tarr)   
3251      dtarr=tarr-0.5       ; convert back to julian day
3252
3253      xmx=max(dtarr,min=xmn)
3254;put a bit of white space around xmx,xmn
3255      dxmxmn=(xmx-xmn)*0.02
3256                xmx=xmx+dxmxmn
3257                xmn=xmn-dxmxmn
3258
3259
3260;date_label = LABEL_DATE(DATE_FORMAT = $
3261;   ['%I:%S', '%H', '%D %M, %Y'])
3262
3263;date_label = LABEL_DATE(DATE_FORMAT = $
3264;   ['%D %M, %Y'])
3265
3266date_label = LABEL_DATE(DATE_FORMAT=['%D-%M','%Y'])
3267
3268
3269
3270   if (stash.txt eq 0) then begin
3271   if (numtimeseries ne 1) then begin                       
3272      plot,dtarr, meants, xstyle=1, linestyle=1, yrange=[ymn,ymx], title=title, $
3273                       ytitle=typestr, xrange=[xmn,xmx], $
3274                       XTICKUNITS=['Time', 'Time'], XTICKFORMAT = ['LABEL_DATE'],YMARGIN=[6,4]
3275;   XTICKFORMAT = ['LABEL_DATE'], $
3276;   XTICKUNITS = ['Day'], $
3277;   XTICKINTERVAL = 4
3278
3279      plot, dtarr, rmsts, xstyle=1, yrange=[ymn,ymx], /noerase, $
3280                       ytitle=typestr, xrange=[xmn,xmx], $
3281                       XTICKUNITS=['Time', 'Time'], XTICKFORMAT = ['LABEL_DATE'],YMARGIN=[6,4]
3282;   XTICKFORMAT = ['LABEL_DATE'], $
3283;   XTICKUNITS = ['Day'], $
3284;   XTICKINTERVAL = 4
3285
3286
3287      xcoord=0.8
3288      ycoord=0.9
3289      ycoord=0.35
3290                ycoord=0.2
3291;     ukmo_legend,xcoord=xcoord,ycoord=ycoord,delta_y=0.05,$
3292;        ['RMS','mean'],linestyle=[0,1],/normal
3293
3294                  ycoord2=ycoord-0.05
3295                  xcoord2=xcoord+0.03
3296                  xcoord3=xcoord+0.05
3297                  plots, [xcoord,xcoord2],[ycoord,ycoord], linestyle=0, /normal
3298                  xyouts, xcoord3, ycoord, 'RMS', /normal
3299                  plots, [xcoord,xcoord2],[ycoord2,ycoord2], linestyle=1, /normal
3300                  xyouts, xcoord3, ycoord2, 'mean',/normal               
3301
3302   endif else begin
3303      plot,dtarr, narr, xstyle=1, title=title, $
3304                       ytitle='number of obs', $
3305                       XTICKUNITS=['Time', 'Time'], XTICKFORMAT = ['LABEL_DATE'],YMARGIN=[6,4]
3306;   XTICKFORMAT = ['LABEL_DATE'], $
3307;   XTICKUNITS = ['Day'], $
3308;   XTICKINTERVAL = 4
3309     
3310   endelse
3311   endif else begin
3312
3313        outfile=stash.outfile+'_timeseries.txt'
3314   print, 'saving data to ',outfile
3315   OPENW, unit, outfile, /get_lun
3316;        printf,unit,'Output timeseries data: '
3317      printf,unit, strtsal
3318                printf,unit, typestr
3319                printf,unit, xrange, yrange
3320                printf,unit, binspday
3321                nel=n_elements(dtarr)
3322                for i=0L,nel-1 do begin
3323                  printf,unit, dtarr(i), narr(i), meants(i), rmsts(i),format='(d18.8,d18.2,d18.8,d18.8)'
3324                endfor
3325        FREE_LUN, unit
3326
3327   endelse
3328
3329;     print,'min/max lonsel: ',min(lonsel),max(lonsel)
3330
3331end
3332
3333PRO timeserieswindow, stash
3334   basepw=WIDGET_BASE(/column)
3335   drawpw = WIDGET_DRAW(basepw, XSIZE=640, YSIZE=480)
3336        buttons=Widget_Base(basepw,row=1)
3337       
3338        tlb = Widget_Base(buttons,Title='Push-Buttons', row=1, Scr_XSize=400,$
3339       /Exclusive)
3340; no_release only sends select events (not release ones)
3341      buttonp1 = Widget_Button(tlb, Value='1 bin per day',uvalue='RADIO1',/no_release)
3342      buttonp2 = Widget_Button(tlb, Value='2 bins',uvalue='RADIO2',/no_release)
3343      buttonp3=  Widget_Button(tlb, Value='4 bins',uvalue='RADIO3',/no_release)
3344      buttonp4 = Widget_Button(tlb, Value='8 bins',uvalue='RADIO4',/no_release)
3345
3346   if (stash.binspday eq 1) then Widget_Control, buttonp1, Set_Button=1
3347   if (stash.binspday eq 2) then Widget_Control, buttonp2, Set_Button=1
3348         if (stash.binspday eq 4) then Widget_Control, buttonp3, Set_Button=1
3349   if (stash.binspday eq 8) then Widget_Control, buttonp4, Set_Button=1
3350 
3351        tlb2 = Widget_Base(buttons,Title='Push-Buttons', row=1, Scr_XSize=100,$
3352       /NonExclusive)
3353   buttonpa1= Widget_Button(tlb2, Value='plot selected period',uvalue='PLOTSEL')
3354
3355
3356   Widget_Control, buttonpa1, Set_Button=stash.plottssel
3357       
3358   buttonpw = WIDGET_BUTTON(buttons, VALUE='Print',uvalue='PRINTTSW')
3359        buttonpw2 = WIDGET_BUTTON(buttons, VALUE='Output',uvalue='OUTPUTTSW')
3360       
3361   WIDGET_CONTROL, basepw, /REALIZE
3362
3363   xyouts,0.2,0.5,'working...',/normal,charsize=2.5
3364   
3365        plottimeseries,stash
3366
3367   WIDGET_CONTROL, basepw, SET_UVALUE=stash 
3368
3369   XMANAGER, 'timeserieswindow', basepw, /NO_BLOCK
3370end
3371
3372PRO timeserieswindow_event, ev
3373  WIDGET_CONTROL, ev.TOP, GET_UVALUE=stash
3374  WIDGET_CONTROL, ev.ID, GET_UVALUE=uval
3375   print,uval
3376   if (uval eq 'PRINTTSW') then begin
3377      ps=1
3378      eps=0
3379      landscape=1
3380      pr2,file='~/timeseries.ps',landscape=landscape,ps=ps,eps=eps,color=1
3381      plottimeseries,stash
3382      prend2,/view
3383   endif
3384        if (uval eq 'OUTPUTTSW') then begin
3385         savetxt=stash.txt
3386                stash.txt=1
3387         plottimeseries,stash
3388                stash.txt=savetxt
3389        endif
3390        if (uval eq 'RADIO1') then begin
3391         stash.binspday=1
3392                plottimeseries,stash         
3393        endif
3394        if (uval eq 'RADIO2') then begin
3395         stash.binspday=2
3396                plottimeseries,stash         
3397        endif
3398        if (uval eq 'RADIO3') then begin
3399         stash.binspday=4
3400                plottimeseries,stash         
3401        endif
3402        if (uval eq 'RADIO4') then begin
3403         stash.binspday=8
3404                plottimeseries,stash         
3405        endif
3406        if (uval eq 'PLOTSEL') then begin
3407         stash.plottssel=1-stash.plottssel
3408                print,'stash.plottssel ',stash.plottssel
3409                WIDGET_CONTROL, ev.TOP, SET_UVALUE=stash
3410                plottimeseries,stash       
3411        endif
3412
3413       
3414end
3415
3416;----------------------------------------
3417; T-S diagram plotting and event handling
3418;----------------------------------------
3419
3420PRO plottsdiagram, stash
3421 
3422      xrange=stash.xrange
3423      yrange=stash.yrange
3424;     dayarr=stash.dayarr
3425      daymin=stash.daymin
3426      daymax=stash.daymax
3427      dayminl=stash.dayminl
3428                daymaxl=stash.daymaxl
3429                xarr=stash.xarr
3430                yarr=stash.yarr
3431                rmdi=stash.rmdi
3432
3433      ; select T and S points to plot
3434      typestr=''
3435            selpoints, stash, lonsel, latsel, innovsel, qcsel, daysel, obstypsel, obnumsel, numsel, typestr, $
3436                  daymin=dayminl, daymax=daymaxl, salinity=1
3437                       
3438                ; select points with salinity
3439               
3440                nel=n_elements(lonsel)
3441               
3442                if (nel gt 0) then begin
3443
3444        bkg=stash.bkg
3445                  bkg2=stash.bkg2
3446                  obs=stash.obs
3447                  obs2=stash.obs2
3448                  dep=stash.dep
3449                  strtarr=0
3450        depmin=stash.depmin
3451        depmax=stash.depmax
3452                 
3453        for i=0L,nel-1 do begin
3454         wh=where(lonsel(i) eq xarr and latsel(i) eq yarr and $
3455                          dep ge depmin and dep le depmax )
3456                        if (wh(0) gt -1) then begin
3457                        if (strtarr eq 0) then begin
3458                          strtarr=1
3459                          bkgsel=bkg(wh)
3460                          bkg2sel=bkg2(wh)
3461                          obssel=obs(wh)
3462                          obs2sel=obs2(wh)
3463                        endif else begin
3464                bkgsel=[bkgsel,bkg(wh)]
3465                          bkg2sel=[bkg2sel,bkg2(wh)]
3466                obssel=[obssel,obs(wh)]
3467                          obs2sel=[obs2sel,obs2(wh)]
3468                        endelse
3469                        endif             
3470                  endfor
3471
3472               ;filter out points with missing data
3473               
3474               wh=where(bkgsel ne rmdi and bkg2sel ne rmdi $
3475                    and obssel ne rmdi and obs2sel ne rmdi)
3476               
3477                if (wh(0) gt -1) then begin
3478                  obssel=obssel(wh)
3479                  bkgsel=bkgsel(wh)
3480                  obs2sel=obs2sel(wh)
3481                  bkg2sel=bkg2sel(wh)             
3482                 
3483                  xmn=min([obssel,bkgsel],max=xmx)
3484                  ymn=min([obs2sel,bkg2sel],max=ymx)
3485               
3486               title='T-S diagram: dep '+string(depmin)+'-'+string(depmax)
3487               
3488      plot,obssel,obs2sel,psym=4, color='FFFFFF'x, xtitle='Temperature',ytitle='Salinity',$
3489                  xrange=[xmn,xmx], yrange=[ymn,ymx],xstyle=1,ystyle=1, title=title
3490                oplot,bkgsel,bkg2sel,psym=5, color='FFCC66'x
3491
3492      xcoord=0.8
3493      ycoord=0.9
3494      ycoord=0.2
3495;     ukmo_legend,xcoord=xcoord,ycoord=ycoord,delta_y=0.05,$
3496;        ['obs','bkg'],psym=[4,5],color=['FFFFFF'x, 'FFCC66'x],/normal
3497
3498                  ycoord2=ycoord-0.05
3499                  xcoord2=xcoord+0.03
3500                  xcoord3=xcoord+0.05
3501                  plots, [xcoord,xcoord2],[ycoord,ycoord], psym=4, color='FFFFFF'x, /normal
3502                  xyouts, xcoord+0.05, ycoord, 'obs', /normal
3503                  plots, [xcoord,xcoord2],[ycoord2,ycoord2], psym=5, color='FFCC66'x, /normal
3504                  xyouts, xcoord+0.05, ycoord2, 'bkg',/normal               
3505
3506
3507               endif
3508               endif
3509
3510end
3511
3512PRO tsdiagramwindow, stash
3513   basepw=WIDGET_BASE(/column)
3514   drawpw = WIDGET_DRAW(basepw, XSIZE=640, YSIZE=480)
3515   buttonpw = WIDGET_BUTTON(basepw, VALUE='Print',uvalue='PRINTTSW')
3516   WIDGET_CONTROL, basepw, /REALIZE
3517   
3518        plottsdiagram,stash
3519
3520   WIDGET_CONTROL, basepw, SET_UVALUE=stash 
3521
3522   XMANAGER, 'tsdiagramwindow', basepw, /NO_BLOCK 
3523end
3524
3525PRO tsdiagramwindow_event, ev
3526  WIDGET_CONTROL, ev.TOP, GET_UVALUE=stash
3527  WIDGET_CONTROL, ev.ID, GET_UVALUE=uval
3528   print,uval
3529   if (uval eq 'PRINTTSW') then begin
3530      ps=1
3531      eps=0
3532      landscape=1
3533      pr2,file='~/tsdiagram.ps',landscape=landscape,ps=ps,eps=eps,color=1
3534      plottsdiagram,stash
3535      prend2,/view
3536   endif
3537end
3538
3539PRO jul_to_dtstr, daysel, datestr, notime=notime
3540
3541;print,'called jul_to_dtstr ',daysel
3542
3543; IDL Julian days start at 12 noon
3544
3545   CALDAT, daysel(0)-0.5, Month, Day, Year, Hour, Minute, Second
3546           datestr=strtrim(fix(year),2)+'/'+ $
3547        strtrim(string(fix(month),format='(i2.2)'),2)+'/'+ $
3548        strtrim(string(fix(day),format='(i2.2)'),2)
3549             
3550             if (n_elements(notime) eq 0) then $
3551               datestr=datestr+' '+$
3552        strtrim(string(fix(hour),format='(i2.2)'),2)+':'+ $
3553        strtrim(string(fix(minute),format='(i2.2)'),2)+':'+ $
3554        strtrim(string(fix(second),format='(i2.2)'),2)
3555
3556;           datedt=jul_to_dt(daysel)
3557;           datedt=datedt(0)
3558;           datestr=strtrim(fix(datedt.year),2)+'/'+ $
3559;       strtrim(string(fix(datedt.month),format='(i2.2)'),2)+'/'+ $
3560;       strtrim(string(fix(datedt.day),format='(i2.2)'),2)+' '+ $
3561;       strtrim(string(fix(datedt.hour),format='(i2.2)'),2)+':'+ $
3562;       strtrim(string(fix(datedt.minute),format='(i2.2)'),2)+':'+ $
3563;       strtrim(string(fix(datedt.second),format='(i2.2)'),2)
3564
3565;print,datestr
3566
3567END
3568
3569;-----------------------------------------------------------------------
3570; MAIN routine
3571; Widget creation routine.
3572;  dataplot, [longitude, latitude, deparr, dayarr, valarr, bkgarr]
3573;+
3574PRO dataplot, indata, rmdi=rmdi, filename=filename, salinity=salinity, batch=batch, $
3575   ps=ps, gif=gif, area=area, typeplot=typeplot, ombtypeplot=ombtypeplot, $
3576        alldays=alldays, depths=depths, $
3577        mx=mx, mn=mn, showmdt=showmdt, obstypeselect=obstypeselect, printobstypes=printobstypes, $
3578        daterange=daterange, timeseries=timeseries, binspday=binspday, plot_bad_obs=plot_bad_obs, $
3579        plot_only_bad_obs=plot_only_bad_obs, $
3580        numtimeseries=numtimeseries, txt=txt, netcdf=netcdf, vertgrad=vertgrad, healthcheck=healthcheck, $
3581        duplicates=duplicates, differences=differences, outfile=outfile, white_on_black=white_on_black, $
3582        bindata=bindata, symscale=symscale, hiresmap=hiresmap, coltable=coltable, $
3583        picsize=picsize, pmulti=pmulti, typeproj=typeproj, dayrange=dayrange, $
3584        notfussy=notfussy, mld=mld, binsize=binsize, filterout=filterout
3585;-
3586   
3587; DJL switch off wave compatibility mode
3588res=execute('waveoff')
3589
3590; string array obstypeselect
3591
3592if (n_elements(obstypeselect) eq 0) then obstypeselect=''
3593if (n_elements(printobstypes) eq 0) then printobstypes=''
3594if (n_elements(filterout) eq 0) then filterout=''
3595
3596if (n_elements(binspday) eq 0) then binspday=4     ; plot timeseries every 6 hours
3597plottssel=0 ; plot selected period of timeseries
3598if (n_elements(plot_bad_obs) eq 0) then plot_bad_obs=0   ; don't plot bad obs
3599
3600if (n_elements(plot_only_bad_obs) eq 0) then plot_only_bad_obs=1  ; plot only bad obs if selected
3601if (n_elements(white_on_black) eq 0) then white_on_black=1
3602
3603if (n_elements(duplicates) eq 0) then duplicates=0 ; plot only duplicates if selected
3604if (n_elements(differences) eq 0) then differences=0  ; plot only differences if selected
3605
3606if (n_elements(numtimeseries) eq 0) then numtimeseries=0 ; plot time series of O-B
3607
3608if (n_elements(txt) eq 0) then txt=0
3609
3610if (n_elements(netcdf) eq 0) then netcdf=0         ; if 1 the output netcdf data
3611
3612print,obstypeselect
3613
3614if (n_elements(salinity) eq 0) then salinity=0
3615density=0
3616if (n_elements(vertgrad) eq 0) then vertgrad=0
3617
3618if (n_elements(mld) eq 0) then mld=0
3619
3620if (n_elements(outfile) eq 0) then outfile="dataplot"
3621
3622if (n_elements(bindata) eq 0) then bindata=0
3623
3624if (n_elements(binsize) eq 0) then binsize=[1.0,1.0]
3625
3626if (n_elements(symscale) eq 0) then symscale=1.0
3627
3628if (n_elements(hiresmap) eq 0) then begin
3629   hires_map=0
3630endif else begin
3631   hires_map=hiresmap
3632endelse
3633
3634if (n_elements(rmdi) eq 0) then rmdi=0
3635
3636if (n_elements(coltable) eq 0) then begin
3637   coltable=-1
3638endif
3639
3640if (n_elements(picsize) ne 2) then picsize=[800,512]  ; default gif resolution
3641
3642if (n_elements(pmulti) eq 0) then pmulti=0      ; default pmulti
3643
3644sz=size(indata)
3645nsz=n_elements(sz)
3646type=sz(nsz-2) ; get type 2 integer, 4 float, 7 string
3647if (type eq 7) then filename=indata
3648
3649if (type ne 7) then begin
3650   nel=n_elements(indata)
3651   nel8=nel/8
3652   numobs=nel8
3653
3654   ; get input values
3655   if (nel gt 0) then begin
3656      indata=reform(indata,numobs,8)
3657      xarr=indata(*,0)
3658      yarr=indata(*,1)
3659      dep=indata(*,2)
3660      dayarrin=indata(*,3)
3661      obs=indata(*,4)
3662      bkg=indata(*,5)
3663      obs2=rmdi
3664      obs3=rmdi
3665      bkg2=rmdi
3666      qcarr=long(indata(*,6))
3667      qcarr2=rmdi
3668      obnum=rmdi
3669      obstypes=long(indata(*,7))
3670      filetype="generic"
3671     
3672      print, 'djl max/min dayarrin ',max(dayarrin), min(dayarrin)
3673      print, 'djl max/min qcarr ',max(qcarr), min(qcarr)
3674     
3675   endif
3676endif else begin
3677   if (n_elements(filename) eq 0) then begin
3678       print,'ERROR: No filename supplied'
3679       print,'call: dataplot, filenamearr'
3680       return
3681   endif
3682endelse
3683
3684;read in data
3685if (n_elements(filename) gt 0) then begin
3686   read_cdfobs, filename, numobs=numobs, Latitudes=yarr, Longitudes=xarr, $
3687         obs=obs, modelvals=bkg, qcs=qcarr, $
3688         ob2=obs2, modelval2=bkg2, qc2=qcarr2, $
3689         ob3=obs3, $
3690         Dates=dayarrin, rmdi=rmdi, $
3691         depths=dep, nodates=0, types=obstypes, $
3692         filetype=filetype, iobs=iobs, jobs=jobs, MDT=MDT, error=error, profilenum=obnum, $
3693         notfussy=notfussy, VarName=varname
3694
3695   if (error eq 1) then begin
3696           numobs=0
3697           dayarrin=0
3698           rmdi=-99999
3699           xarr=rmdi
3700           yarr=rmdi
3701           obs=rmdi
3702           bkg=rmdi
3703           qcarr=1
3704           dayarrin=rmdi
3705           dayarr=rmdi
3706           dep=rmdi
3707           filetype=""
3708           obstypes=rmdi
3709           obnum=0
3710   endif
3711
3712   print,'error ',error, numobs
3713
3714   if (numobs gt 0) then begin
3715        print, 'numobs ',numobs
3716
3717;setup obnum if it is not produced by read_cdfobs
3718
3719   if (n_elements(obnum) ne numobs and numobs gt 0) then begin
3720      obnum=lindgen(n_elements(numobs))   ; generate an array of observation numbers
3721                              ; starting with zero
3722        endif
3723
3724
3725
3726;if gt 180 then make longitudes negative
3727
3728     wh=where(xarr ge 180 and xarr le 360)
3729     if (wh(0) gt -1) then begin
3730       xarr(wh)=xarr(wh)-360
3731     endif
3732
3733
3734;calculate vertical gradient
3735;
3736;    if (filetype eq "Prof") then $
3737;    calcvertgrad, xarr, yarr, dayarrin, obs, bkg, qcarr, obs2,bkg2, qcarr2, dep, rmdi
3738
3739; adjust qc values for profile obs data (1 = good)
3740     wh=where(bkg ne rmdi)
3741          if (n_elements(bkg2) gt 0) then begin
3742            wh2=where(bkg2 ne rmdi)
3743       if (filetype eq "Prof" and wh(0) eq -1 and wh2(0) eq -1) then begin
3744              qcarr=qcarr-1
3745              qcarr2=qcarr2-1
3746       endif
3747     endif
3748 
3749          if (keyword_set(showmdt) eq 1) then begin
3750       bkg=MDT
3751          endif
3752
3753
3754     print,'rmdi ',rmdi
3755
3756     print,'file read in'
3757
3758spawn,'echo part 1 `date`'
3759
3760     if (n_elements(rmdi) eq 0) then rmdi=-1.07374e+09
3761
3762     qcarr=fix(qcarr)
3763
3764   endif    ; numobs
3765
3766   print,'endif numobs'
3767
3768; select area
3769
3770        sz=size(area)
3771        nsz=n_elements(sz)
3772        type=sz(nsz-2)  ; get type 2 integer, 4 float, 7 string
3773        if (type eq 7) then areasel, area, xrange, yrange, success
3774        if (type eq 4 or type eq 2) then begin
3775     if (nsz eq 4) then begin
3776             xrange=[area(0),area(2)]
3777             yrange=[area(1),area(3)]
3778          endif
3779   endif       
3780
3781; set missing obs types to missing data
3782        if (n_elements(obs2) eq 0) then begin
3783          obs2=rmdi
3784          bkg2=rmdi
3785          qcarr2=rmdi
3786       endif
3787        if (n_elements(obs3) eq 0) then begin
3788          obs3=rmdi
3789          bkg3=rmdi
3790          qcarr3=rmdi
3791       endif
3792
3793endif    ; filename
3794
3795if (n_elements(varname) eq 0) then varname=''
3796
3797     if (n_elements(dayarrin) ne numobs) then dayarrin=replicate(1,numobs)
3798;      dep=replicate(0,numobs)
3799
3800
3801     if (size(dayarrin,/type) eq 8) then begin
3802       dayarr=dayarrin.julian
3803          endif else begin
3804       dayarr=dayarrin
3805            dayarr=dayarr+0.5       ; dataplot prefers each day to be a different integer
3806          endelse
3807
3808print,'xxx'
3809
3810;set contour range mn, mx
3811if (n_elements(mx) eq 0) then fmx=rmdi else fmx=mx
3812if (n_elements(mn) eq 0) then fmn=rmdi else fmn=mn
3813mx=fmx
3814mn=fmn
3815
3816
3817
3818  ; Define a monochrome image array for use in the application.
3819;  READ_PNG, FILEPATH('mineral.png', $
3820;    SUBDIR=['examples', 'data']), image
3821 
3822  ; Place the image array in a pointer heap variable, so we can
3823  ; pass the pointer to the event routine rather than passing the
3824  ; entire image array.
3825;  imagePtr=PTR_NEW(image, /NO_COPY)
3826 
3827  ; Retrieve the size information from the image array.
3828;  im_size=SIZE(*imagePtr)
3829
3830;---------------
3831;Setup defaults
3832;---------------
3833
3834  busy=0
3835  wh=where(dep ge 0)
3836  depminl=0.
3837  depmaxl=1.
3838  if( wh(0) gt -1) then depminl=min(dep(wh),max=depmaxl)
3839  depscl=fix((depmaxl-depminl)/100.)
3840  if (depscl lt 1) then depscl=1
3841  depmin=depminl
3842;  depmax=depmin+depscl 
3843  depmax=depmaxl
3844  pmultinum=0
3845
3846if (n_elements(depths) eq 2) then begin
3847
3848 depmin=max([depminl,depths(0)])
3849 depmax=min([depmaxl,depths(1)])
3850endif
3851
3852   print,'n_elements(typeplot) ',n_elements(typeplot)
3853
3854   if (n_elements(typeplot) eq 0) then typeplot=1
3855;   if (n_elements(where(bkg ne bkg(0))) gt 1) then begin
3856;     typeplot=3
3857;   endif else begin
3858;     typeplot=4
3859;   endelse
3860;   endif
3861   if (typeplot gt 4) then typeplot=4
3862   if (typeplot lt 1) then typeplot=1
3863   
3864   if (n_elements(ombtypeplot) eq 0) then begin
3865    if (n_elements(where(bkg ne bkg(0))) gt 1) then begin
3866      ombtypeplot=1
3867    endif else begin
3868      ombtypeplot=2
3869    endelse
3870   endif
3871   if (ombtypeplot gt 3) then ombtypeplot=3
3872   if (ombtypeplot lt 1) then ombtypeplot=1
3873
3874      print,'typeplot ',typeplot, ' ombtypeplot ',ombtypeplot
3875   
3876    if (n_elements(typeproj) eq 0) then typeproj=1
3877
3878   xrangedef=[-110.,40.]
3879   yrangedef=[-50.,80.]
3880   wh=where(xarr lt xrangedef(0) or xarr gt xrangedef(1) or $
3881      yarr lt yrangedef(0) or yarr gt yrangedef(1),count)
3882   if (count gt 0) then begin
3883      xrangedef=[-180.,180.]
3884      yrangedef=[-90.,90.]
3885   endif
3886
3887   if (n_elements(xrange) eq 0) then begin
3888   xrange=xrangedef
3889   yrange=yrangedef
3890        endif
3891   symsize=1.0
3892
3893;  map_file=""
3894
3895;  plot_bad_obs=0
3896
3897  dayshi=-10.d/86400.      ; small shift to group 0:00 hours with the previous day
3898;  dayshi=-0.1d
3899
3900  wh=where(dayarr gt 0 and dayarr lt 9000000)
3901;  print,dayarr
3902;  print,wh(0)
3903  if (wh(0) gt -1) then begin
3904    daymin=min(dayarr(wh)+dayshi,max=daymax)
3905  endif else begin
3906    daymax=0
3907    daymin=0
3908  endelse
3909;  daymax=maxday
3910;  daymin=minday
3911  dayminl=daymin
3912  daymaxl=daymax
3913 
3914
3915  print,'** dayminl daymaxl ', dayminl, daymaxl
3916
3917; ** Health check file   
3918   if (keyword_set(healthcheck)) then begin
3919      OPENW, unit, outfile+'_health.txt', /get_lun
3920         printf,unit,' Health check data '
3921                printf,unit,'-------------------'
3922                printf,unit,' Number of files ',n_elements(filename)
3923                printf,unit,' List of files '
3924      printf,unit,filename
3925                printf,unit,' Filetype ',filetype
3926                printf,unit,' Number of observations ',numobs
3927                printf,unit,' Date range in julian days ',dayminl, daymaxl,format='(a,2f20.8)'
3928      jul_to_dtstr,dayminl,datestr1
3929                jul_to_dtstr,daymaxl,datestr2
3930      printf,unit,' Date range ',datestr1,' - ',datestr2
3931                mxo=0 & mno=0
3932                mxb=0 & mnb=0
3933                mxob=0 & mnob=0
3934                rmsomb=0 & meanomb=0
3935;                info,obs
3936;                print,obs
3937;                info,bkg
3938;                print,bkg
3939;                print,'rmdi: ',rmdi
3940;                wh=where(obs ne rmdi and qcarr lt 1,counto)
3941; check for missing data within a tolerance
3942      wh=where(abs((obs-rmdi)/rmdi) gt 0.01 and qcarr lt 1, counto)
3943;                print,'min ',min(abs((bkg(wh)-rmdi)/rmdi))
3944;                print,'min ',min(bkg(wh)), min(obs(wh))
3945;                print,'min ',max(bkg(wh)), max(obs(wh))
3946                if (counto gt 0) then mxo=max(obs(wh),min=mno)
3947;               wh=where(bkg ne rmdi and qcarr lt 1,countb)
3948      wh=where(abs((bkg-rmdi)/rmdi) gt 0.01 and qcarr lt 1, countb)
3949                if (countb gt 0) then mxb=max(bkg(wh),min=mnb)
3950;                wh=where(obs ne rmdi and bkg ne rmdi and qcarr lt 1,countob)
3951      wh=where(abs((obs-rmdi)/rmdi) gt 0.01 and $
3952                  abs((bkg-rmdi)/rmdi) gt 0.01 and qcarr lt 1, countob)
3953                if (countob gt 0) then mxob=max(obs(wh)-bkg(wh),min=mnob)
3954; calculate rms / mean
3955      if (countob gt 0) then begin
3956         x=total(obs(wh)-bkg(wh))
3957                   x2=total((obs(wh)-bkg(wh))^2)
3958                   nel=n_elements(wh)
3959                   meanomb=x/nel
3960                   rmsomb=sqrt(x2/nel)
3961                endif
3962 
3963                printf,unit,' Max/min obs ',mxo,mno
3964                printf,unit,' Max/min bkg ',mxb,mnb
3965                printf,unit,' Max/min obs-bkg ',mxob,mnob
3966                printf,unit,' RMS/mean obs-bkg ',rmsomb,meanomb
3967                printf,unit,' Number of good observations ',counto
3968                printf,unit,' Number of bad observations ', numobs-counto
3969
3970                mxo=0 & mno=0
3971                mxb=0 & mnb=0
3972                mxob=0 & mnob=0
3973                rmsomb=0 & meanomb=0
3974                wh=where(obs2 ne rmdi and qcarr2 lt 1,counto)
3975                if (counto gt 0) then mxo=max(obs2(wh),min=mno)
3976                wh=where(bkg2 ne rmdi and qcarr2 lt 1,countb)
3977                if (countb gt 0) then mxb=max(bkg2(wh),min=mnb)
3978                wh=where(obs2 ne rmdi and bkg2 ne rmdi and qcarr2 lt 1,countob)
3979                if (countob gt 0) then mxob=max(obs2(wh)-bkg2(wh),min=mnob)
3980; calculate rms / mean
3981      if (countob gt 0) then begin
3982         x=total(obs2(wh)-bkg2(wh))
3983                   x2=total((obs2(wh)-bkg2(wh))^2)
3984                   nel=n_elements(wh)
3985                   meanomb=x/nel
3986                   rmsomb=sqrt(x2/nel)
3987                endif
3988               
3989                printf,unit,' Max/min obs2 ',mxo,mno
3990                printf,unit,' Max/min bkg2 ',mxb,mnb
3991                printf,unit,' Max/min obs2-bkg2 ',mxob,mnob
3992                printf,unit,' RMS/mean obs2-bkg2 ',rmsomb,meanomb               
3993                printf,unit,' Number of good observations2 ',counto
3994                printf,unit,' Number of bad observations2 ', numobs-counto
3995
3996           FREE_LUN,unit
3997        endif
3998
3999        daymax=daymaxl        ; default to show everything
4000        if (keyword_set(alldays)) then daymax=daymaxl         ;print alldays
4001        if (n_elements(daterange) eq 2) then begin
4002         daymin=daterange(0)
4003                if (daterange(0) lt dayminl) then daymin=dayminl
4004                if (daterange(0) gt daymaxl) then daymin=daymaxl
4005                if (daterange(1) gt dayminl and daterange(1) le daymaxl) then begin
4006                  daymax=daterange(1)
4007                endif else begin
4008                daymax=daymaxl               
4009                endelse
4010        endif
4011        if (n_elements(daterange) eq 1) then begin    ; select number of days before end or after beginning
4012               print,'daymin b4: ',daymin, daymax
4013                daymax=daymaxl
4014                daymin=dayminl
4015                if (daterange(0) lt 0) then begin
4016                  if (daterange(0) eq -9999) then begin
4017                           daymin=daymax
4018                        endif else begin
4019                     daymin=daymax+daterange(0)
4020                        endelse
4021                endif
4022                if (daterange(0) gt 0) then daymax=daymin+daterange(0)
4023        endif                                 
4024
4025; select day from the start day
4026   if (n_elements(dayrange) eq 1) then begin
4027         daymin=dayminl+dayrange
4028                daymax=daymin
4029        endif
4030
4031   if (n_elements(dayrange) eq 2) then begin
4032         daymin=dayminl+dayrange(0)
4033                daymax=dayminl+dayrange(1)
4034        endif
4035
4036
4037; prevent error where only one day plotted 
4038   daymaxsl=daymaxl
4039        dayminsl=dayminl
4040   if (daymaxsl le dayminsl+1) then daymaxsl=dayminsl+1
4041; stop day going out of range
4042   if (daymax gt daymaxl) then daymax=daymaxl
4043        if (daymax lt dayminl) then daymax=dayminl
4044        if (daymin lt dayminl) then daymin=dayminl
4045        if (daymin gt daymaxl) then daymin=daymaxl
4046
4047   print,'daymin af: ',daymin, daymax
4048
4049  sym=1
4050   
4051;-------------
4052;Setup window
4053;-------------
4054
4055if (not keyword_set(batch)) then begin
4056  ; Create a base widget to hold the application.
4057;  base0=WIDGET_BASE()
4058
4059  base = WIDGET_BASE(/COLUMN,mbar=bar)
4060;   base = WIDGET_BASE(/row)
4061
4062 
4063  ; setup menu bar
4064 
4065  menu1 = WIDGET_BUTTON(bar, VALUE='Areas', /MENU)
4066button1 = WIDGET_BUTTON(menu1, VALUE='Global', uvalue='Global')
4067button2 = WIDGET_BUTTON(menu1, VALUE='Arctic', uvalue='Arctic')
4068button2a = WIDGET_BUTTON(menu1, VALUE='Atlantic', uvalue='Atlantic')
4069button3 = WIDGET_BUTTON(menu1, VALUE='N Atl', uvalue='N Atl')
4070button4 = WIDGET_BUTTON(menu1, VALUE='Trop Atl', uvalue='Trop Atl')
4071button5 = WIDGET_BUTTON(menu1, VALUE='S Atl', uvalue='S Atl')
4072button5a = WIDGET_BUTTON(menu1, VALUE='Pacific', uvalue='Pacific')
4073button6 = WIDGET_BUTTON(menu1, VALUE='N Pac', uvalue='N Pac')
4074button7 = WIDGET_BUTTON(menu1, VALUE='Trop Pac', uvalue='Trop Pac')
4075button8 = WIDGET_BUTTON(menu1, VALUE='S Pac', uvalue='S Pac')
4076button9 = WIDGET_BUTTON(menu1, VALUE='Indian', uvalue='Indian')
4077button10 = WIDGET_BUTTON(menu1, VALUE='S Ocean', uvalue='S Ocean')
4078button11 = WIDGET_BUTTON(menu1, VALUE='Med', uvalue='Med')
4079button12 = WIDGET_BUTTON(menu1, VALUE='N West Shelf', uvalue='NWS')
4080button13 = WIDGET_BUTTON(menu1, VALUE='Input Area', uvalue='Input Area')
4081
4082  menu2 = WIDGET_BUTTON(bar, VALUE='Plot', /MENU)
4083button1 = WIDGET_BUTTON(menu2, VALUE='Timeseries', uvalue='Timeseries')
4084button1a =  WIDGET_BUTTON(menu2, VALUE='Num timeseries', uvalue='Num timeseries')
4085button2 = WIDGET_BUTTON(menu2, VALUE='T-S diagram', uvalue='TS diagram')
4086
4087  menu1a = WIDGET_BUTTON(bar, VALUE='Find',/MENU)
4088button1a= WIDGET_BUTTON(menu1a, VALUE='Worst points', uvalue='Worst points')
4089button1a1 = WIDGET_BUTTON(menu1a, VALUE='Filter type', uvalue='Filter type')
4090
4091  menu3 = WIDGET_BUTTON(bar, VALUE='Config', /MENU)
4092loresmap = WIDGET_BUTTON(menu3, VALUE='Low res map', uvalue='Low res map', /checked_menu)
4093hiresmap = WIDGET_BUTTON(menu3, VALUE='Hi res map', uvalue='Hi res map', /checked_menu)
4094
4095button3 = WIDGET_BUTTON(menu3, VALUE='Square psym', uvalue='Square psym')
4096button4 = WIDGET_BUTTON(menu3, VALUE='Round psym', uvalue='Round psym')
4097
4098pltonbado = WIDGET_BUTTON(menu3, VALUE='Plot only bad obs', uvalue='Plot only bad obs', /checked_menu)
4099
4100whtonblack = WIDGET_BUTTON(menu3, VALUE='White on black', uvalue='White on black', /checked_menu)
4101
4102incsym = WIDGET_BUTTON(menu3, VALUE='Psym size+ ', uvalue='incsym')
4103decsym = WIDGET_BUTTON(menu3, VALUE='Psym size- ', uvalue='decsym')
4104resetsym = WIDGET_BUTTON(menu3, VALUE='Psym size reset ', uvalue='resetsym')
4105
4106vertgradmenu = WIDGET_BUTTON(menu3, VALUE='Vert gradient', uvalue='vertgrad', /checked_menu)
4107button12 = WIDGET_BUTTON(menu3, VALUE='Input max/min', uvalue='input max/min')
4108
4109
4110  menu4 = WIDGET_BUTTON(bar, VALUE='Help', /MENU)
4111help1 = WIDGET_BUTTON(menu4, VALUE='Info', uvalue='Info')
4112
4113if (n_elements(filename) gt 0) then begin
4114  menu5 = WIDGET_BUTTON(bar, VALUE=filename[0])
4115nel=n_elements(filename)
4116nel2=nel
4117nelmx=32
4118if (nel2 gt nelmx) then nel2=nelmx
4119for i=0L, nel2-1 do begin
4120  id = WIDGET_BUTTON(menu5, VALUE=filename[i])
4121endfor
4122if (nel gt nelmx) then id=WIDGET_BUTTON(menu5, VALUE='... etc ...')
4123endif
4124
4125Widget_Control, loresmap, Set_Button=1
4126Widget_Control, vertgradmenu, Set_Button=vertgrad
4127
4128Widget_Control, pltonbado, Set_Button=plot_only_bad_obs
4129Widget_Control, whtonblack, Set_Button=white_on_black
4130 
4131  ; Create a draw widget based on the size of the image, and
4132  ; set the MOTION_EVENTS keyword so that events are generated
4133  ; as the cursor moves across the image. Setting the BUTTON_EVENTS
4134  ; keyword rather than MOTION_EVENTS would require the user to
4135  ; click on the image before an event is generated.
4136
4137;  im_size=[2,1024,800]
4138;  im_size=[2,800,800]
4139;  im_size=[2,800,720]
4140;  im_size=[2,1024,720]
4141        im_size=[2,1024,680]
4142;  im_size=[2,800,640]
4143
4144  draw = WIDGET_DRAW(base, XSIZE=im_size[1], YSIZE=im_size[2], $
4145    /BUTTON_EVENTS,/WHEEL_EVENTS, uvalue='MAPWINDOW')
4146
4147;  base2=widget_base(GROUP_LEADER=base,/column)
4148
4149;  droplist = WIDGET_COMBOBOX(base2,value=string([indgen(101)]),$
4150;     uvalue='LEVELLIST', FRAME=2)
4151
4152;  button = WIDGET_BUTTON(base2, value='test',uvalue='test')
4153
4154
4155;  print,'depminl/maxl ',depminl, depmaxl
4156;        print,'value depmin ',depmin
4157;        print,'value depmax ',depmax
4158
4159depmaxlsl=depmaxl
4160depminlsl=depminl
4161if (depmaxlsl eq depminlsl) then depmaxlsl=depmaxlsl+1
4162print, 'depminlsl, depmaxlsl ',depminlsl, depmaxlsl
4163  slider = WIDGET_SLIDER(base,maximum=depmaxlsl,minimum=depminlsl,scroll=depscl,$
4164   value=depmin,uvalue='LEVELCHOICE', ysize=32)
4165
4166  sliderb = WIDGET_SLIDER(base,maximum=depmaxlsl,minimum=depminlsl, scroll=depscl,$
4167        value=depmax,uvalue='LEVELCHOICEB', ysize=32)
4168
4169
4170  jul_to_dtstr,daymin,dayminstr,/notime
4171  slider1label = WIDGET_LABEL(base, VALUE="min date: "+dayminstr, xsize=800,ysize=14)
4172  slider1 = WIDGET_SLIDER(base,minimum=dayminsl,maximum=daymaxsl,scroll=1,$
4173   value=daymin,uvalue='DATERANGE1',/suppress_value,/drag)
4174
4175  jul_to_dtstr,daymax,daymaxstr,/notime
4176  slider2label = WIDGET_LABEL(base, VALUE="max date: "+daymaxstr, xsize=800,ysize=14)
4177  slider2 = WIDGET_SLIDER(base,minimum=dayminsl,maximum=daymaxsl,scroll=1,$
4178   value=daymax,uvalue='DATERANGE2',/suppress_value,/drag)
4179
4180
4181
4182;  slider1 = WIDGET_SLIDER(base,minimum=dayminl,maximum=daymaxl,scroll=1,$
4183;     value=daymin,uvalue='DATERANGE1',PRO_SET_VALUE="datesliderset",/drag)
4184;
4185;  slider2 = WIDGET_SLIDER(base,minimum=dayminl,maximum=daymaxl,scroll=1,$
4186;     value=daymax,uvalue='DATERANGE2',PRO_SET_VALUE="datesliderset",/drag)
4187   
4188  ; Labels for stats
4189   
4190  labt1='Lon: '
4191  labt2='Lat: '
4192  labt3='Value: '
4193  labt4='QC: '
4194  labt5='Date: '
4195  labt6='Value: '
4196  labt7=' '
4197  labt8='Type: '
4198  labt9='obnum: '
4199  labt10='visible date range: '
4200 
4201  ; Create label widgets to hold the cursor position and
4202  ; Hexadecimal value of the pixel under the cursor.
4203  labelb=Widget_Base(base,row=1)
4204 
4205  labsiz=105
4206  labsizl1=150
4207  labsizlg=150
4208  labsizl2=135
4209  label7=0
4210  lb_ysize=14
4211  label1 = WIDGET_LABEL(labelb, XSIZE=labsiz, $
4212    VALUE=labt1, ysize=lb_ysize)
4213  label2 = WIDGET_LABEL(labelb, XSIZE=labsiz, $
4214    VALUE=labt2, ysize=lb_ysize)
4215  label3 = WIDGET_LABEL(labelb, XSIZE=labsizl1, $
4216    VALUE=labt3, ysize=lb_ysize)
4217  label4 = WIDGET_LABEL(labelb, XSIZE=labsiz, $
4218    VALUE=labt4, ysize=lb_ysize)
4219  label5 = WIDGET_LABEL(labelb, XSIZE=labsizlg, $
4220    VALUE=labt5, ysize=lb_ysize)
4221  label6 = WIDGET_LABEL(labelb, XSIZE=labsiz, $
4222    VALUE=labt6, ysize=lb_ysize)
4223;  label7 = WIDGET_LABEL(labelb, XSIZE=labsiz, $
4224;    VALUE=labt7, ysize=lb_ysize)
4225  label8 = WIDGET_LABEL(labelb, XSIZE=labsizl2, $
4226    VALUE=labt8, ysize=lb_ysize)
4227  label9 = WIDGET_LABEL(labelb, XSIZE=labsiz, $
4228    VALUE=labt9, ysize=lb_ysize)
4229
4230  labelb2 = Widget_Base(base,row=1)
4231  label10 = WIDGET_LABEL(labelb2, XSIZE=1024, $
4232    VALUE=labt10, ysize=lb_ysize)
4233
4234   tlba = Widget_Base(base,row=1)
4235
4236   t1b = Widget_Base(tlba,Title='Push-Buttons', row=1, Scr_XSize=200,$
4237   /Exclusive)
4238   button1 = Widget_Button(t1b, Value='mean',uvalue='RADIO1',/no_release)
4239   button2 = Widget_Button(t1b, Value='rms',uvalue='RADIO2',/no_release)
4240   button3 = Widget_Button(t1b, Value='sd',uvalue='RADIO3',/no_release)
4241   button4 = Widget_Button(t1b, Value='ms',uvalue='RADIO4',/no_release)
4242
4243   if (typeplot eq 1) then Widget_Control, button1, Set_Button=typeplot
4244   if (typeplot eq 2) then Widget_Control, button2, Set_Button=typeplot
4245   if (typeplot eq 3) then Widget_Control, button3, Set_Button=typeplot
4246   if (typeplot eq 4) then Widget_Control, button4, Set_Button=typeplot
4247 
4248   t1c = Widget_Base(tlba,Title='Push-Buttons', row=1, Scr_XSize=200,$
4249   /Exclusive)
4250   button5 = Widget_Button(t1c, Value='obs - bkg', uvalue='RADIO5',/no_release)
4251   button6 = Widget_Button(t1c, Value='obs', uvalue='RADIO6',/no_release)
4252   button7 = Widget_Button(t1c, Value='bkg', uvalue='RADIO7',/no_release)
4253
4254   if (ombtypeplot eq 1) then Widget_Control, button5, /Set_Button
4255   if (ombtypeplot eq 2) then Widget_Control, button6, /Set_Button
4256   if (ombtypeplot eq 3) then Widget_Control, button7, /Set_Button
4257
4258;   tlb = Widget_Base(tlba,Title='Push-Buttons', row=1, Scr_XSize=400,$
4259;   /Exclusive)
4260;; no_release only sends select events (not release ones)
4261;   button1 = Widget_Button(tlb, Value='rms Obs - Bkg',uvalue='RADIO1',/no_release)
4262;   button2 = Widget_Button(tlb, Value='sum(Obs - Bkg)^2',uvalue='RADIO2',/no_release)
4263;   button3=  Widget_Button(tlb, Value='Obs - Bkg',uvalue='RADIO3',/no_release)
4264;   button4 = Widget_Button(tlb, Value='Obs',uvalue='RADIO4',/no_release)
4265;   button5 = Widget_Button(tlb, Value='Bkg',uvalue='RADIO5',/no_release)
4266;
4267;   if (typeplot eq 1) then Widget_Control, button1, Set_Button=typeplot
4268;   if (typeplot eq 2) then Widget_Control, button2, Set_Button=typeplot
4269;   if (typeplot eq 3) then Widget_Control, button3, Set_Button=typeplot
4270;   if (typeplot eq 4) then Widget_Control, button4, Set_Button=typeplot
4271;   if (typeplot eq 5) then Widget_Control, button5, Set_Button=typeplot
4272 
4273   tlb1 = Widget_Base(tlba,Title='Push-Buttons', row=1, Scr_XSize=220,$
4274   /Exclusive)
4275   button01 = Widget_Button(tlb1, Value='Lat/Lon',uvalue='RADIO01',/no_release)
4276   button02 = Widget_Button(tlb1, Value='N Polar',uvalue='RADIO02',/no_release)
4277   button03 = Widget_Button(tlb1, Value='S Polar',uvalue='RADIO03',/no_release)
4278
4279    Widget_Control, button01, Set_Button=typeproj
4280
4281
4282   tlb1a = Widget_Base(tlba,Title='Push-Buttons', row=1, Scr_XSize=250,$
4283   /NonExclusive)
4284   button001 = Widget_Button(tlb1a, Value='Plot Bad Obs', uvalue='RADIO001')
4285
4286   Widget_Control, button001, Set_Button=plot_bad_obs
4287
4288
4289   if filetype eq 'CRT' then begin
4290
4291      button002 = Widget_Button(tlb1a, Value='Meridional', uvalue='RADIOSAL')
4292      Widget_Control, button002, Set_Button=salinity
4293
4294      button003 = Widget_Button(tlb1a, Value='Speed', uvalue='RADIODENSITY')
4295      Widget_Control, button003, Set_Button=density
4296
4297   endif else begin
4298;         tlb1b = Widget_Base(tlba,Title='Push-Buttons', row=1, Scr_XSize=90,$
4299;         /NonExclusive)
4300      button002 = Widget_Button(tlb1a, Value='Salinity', uvalue='RADIOSAL')
4301      Widget_Control, button002, Set_Button=salinity
4302;         tlb1c = Widget_Base(tlba,Title='Push-Buttons', row=1, Scr_XSize=90,$
4303;          /NonExclusive)
4304      button003 = Widget_Button(tlb1a, Value='Density', uvalue='RADIODENSITY')
4305      Widget_Control, button003, Set_Button=density
4306   endelse
4307
4308 
4309   tlb2= Widget_Base(tlba,row=1)
4310  button = WIDGET_BUTTON(tlb2, VALUE='Print',uvalue="PRINT")
4311  button = WIDGET_BUTTON(tlb2, VALUE='Save',uvalue="SAVE")
4312   
4313  button = WIDGET_BUTTON(tlb2, VALUE='Done',uvalue="DONE")
4314
4315
4316
4317  ; Realize the widget hierarchy.
4318  WIDGET_CONTROL, base, /REALIZE
4319;  WIDGET_CONTROL, base2, /REALIZE
4320;   WIDGET_CONTROL, tlb, /REALIZE
4321
4322  ; Retrieve the widget ID of the draw widget. Note that the widget
4323  ; hierarchy must be realized before you can retrieve this value.
4324  WIDGET_CONTROL, draw, GET_VALUE=drawID
4325
4326 
4327  ; Create an anonymous array to hold the image data and widget IDs
4328  ; of the label widgets.
4329;  stash = { imagePtr:imagePtr, label1:label1, label2:label2, $
4330;            label3:label3, $
4331;      xarr:xarr, yarr:yarr, dep:dep, val:val }
4332
4333
4334spawn,'echo part 2 `date`'
4335
4336
4337  stash = { label1:label1, label2:label2, $
4338            label3:label3, label4:label4, $
4339            label5:label5, label6:label6, $
4340            label7:label7, label8:label8, $
4341            label9:label9, label10:label10, $
4342            labt1:labt1, labt2:labt2, $
4343            labt3:labt3, labt4:labt4, $
4344            labt5:labt5, labt6:labt6, $
4345            labt7:labt7, labt8:labt8, $
4346            labt9:labt9, labt10:labt10, $
4347            slider1label:slider1label, $
4348            slider2label:slider2label, $
4349            drawID:drawID, draw:draw, base:base, $
4350            im_size:im_size, pixID:0, $
4351            xcorn:0, ycorn:0, xdatacorn:0.0, ydatacorn:0.0, $
4352            slider1:slider1, slider2:slider2, $
4353            slider:slider, sliderb:sliderb, $
4354            hiresmap:hiresmap, loresmap:loresmap, vertgradmenu:vertgradmenu, $
4355            pltonbado:pltonbado, plot_only_bad_obs:plot_only_bad_obs, $
4356            whtonblack:whtonblack, white_on_black:white_on_black, $
4357            duplicates:duplicates, differences:differences, $
4358            bindata:bindata, binsize:binsize, $
4359            filetype:filetype, $
4360            xarr:xarr, yarr:yarr, dep:dep, dayarr:dayarr, $
4361            obs:obs, bkg:bkg, qcarr:qcarr, obstypes: obstypes, $
4362            obs2:obs2, obs3:obs3, bkg2:bkg2, qcarr2:qcarr2, $
4363            obnum:obnum, $
4364            depmin:depmin, typeplot:typeplot, ombtypeplot:ombtypeplot, $
4365            typeproj:typeproj, $
4366            hires_map:hires_map, $
4367            plot_bad_obs:plot_bad_obs, salinity: salinity, $
4368            density:density, vertgrad:vertgrad, mld:mld, $
4369            daymin:daymin, daymax:daymax, dayminl:dayminl, daymaxl:daymaxl, $
4370            dayshi:dayshi, $
4371            depminl:depminl, depmaxl:depmaxl, depmax:depmax, depscl:depscl, $
4372            xrange:xrange, yrange:yrange, $
4373            fmx:fmx, fmn:fmn, mx:mx, mn:mn, $
4374            symsize:symsize, sym:sym, picsize:picsize, $
4375            pmulti:pmulti, pmultinum:pmultinum, $
4376            outfile:outfile, $
4377            xrangedef:xrangedef, yrangedef:yrangedef, rmdi:rmdi, $
4378            obstypeselect:obstypeselect, printobstypes:printobstypes, $
4379            plottssel:plottssel, $
4380            binspday:binspday, numtimeseries:numtimeseries, txt:txt, netcdf:netcdf, $
4381            busy:busy, symscale:symscale, coltable:coltable, $
4382            filterout:filterout, varname:varname }
4383       
4384  ; Set the user value of the top-level base widget equal to the
4385  ; 'stash' array.
4386  WIDGET_CONTROL, base, SET_UVALUE=stash
4387 
4388  ; Make the draw widget the current IDL drawable area.
4389  WSET, drawID
4390
4391spawn,'echo part 3 `date`'
4392
4393  ; Draw the image into the draw widget.
4394  plotpoints, stash
4395
4396; plotpoints updates mx/mn values make sure these are included
4397  WIDGET_CONTROL, base, SET_UVALUE=stash
4398
4399spawn,'echo part 4 `date`'
4400
4401  ; Call XMANAGER to manage the widgets.
4402  XMANAGER, 'dataplot', base, /NO_BLOCK
4403;  XMANAGER, 'dataplot', base
4404endif else begin
4405
4406;do batch
4407
4408print,"batch"
4409
4410drawID=-1
4411
4412print, 'filetype: ',filetype
4413
4414  stash = { drawID:drawID,$       
4415            filetype:filetype, $
4416            xarr:xarr, yarr:yarr, dep:dep, dayarr:dayarr, $
4417            obs:obs, bkg:bkg, qcarr:qcarr, obstypes: obstypes, $
4418            obs2:obs2, bkg2:bkg2, qcarr2:qcarr2, $
4419            obnum:obnum, $
4420            depmin:depmin, typeplot:typeplot, ombtypeplot:ombtypeplot, $
4421            typeproj:typeproj, $
4422            hires_map:hires_map, $
4423            plot_only_bad_obs:plot_only_bad_obs, $
4424            white_on_black:white_on_black, $
4425            duplicates:duplicates, differences:differences, $
4426            bindata: bindata, binsize:binsize, $
4427            plot_bad_obs:plot_bad_obs, salinity: salinity, $
4428            density:density, vertgrad:vertgrad, mld:mld, $
4429            daymin:daymin, daymax:daymax, dayminl:dayminl, daymaxl:daymaxl, $
4430            dayshi:dayshi, $           
4431            depminl:depminl, depmaxl:depmaxl, depmax:depmax, depscl:depscl, $
4432            xrange:xrange, yrange:yrange, $
4433            fmx:fmx, fmn:fmn, mx:mx, mn:mn, $
4434            symsize:symsize, sym:sym, picsize:picsize, $
4435            pmulti:pmulti, pmultinum:pmultinum, $
4436            outfile:outfile, $
4437            xrangedef:xrangedef, yrangedef:yrangedef, rmdi:rmdi, $
4438            obstypeselect:obstypeselect, printobstypes:printobstypes, $
4439            plottssel:plottssel, $
4440            binspday:binspday, numtimeseries:numtimeseries, txt:txt, netcdf:netcdf, $
4441            busy:busy, symscale:symscale, coltable:coltable, $
4442            filterout:filterout, varname:varname }
4443
4444if (keyword_set(ps)) then begin
4445   ps=1
4446   eps=0
4447   landscape=1
4448   pr2,file=stash.outfile+'.ps',landscape=landscape,ps=ps,eps=eps,color=1 
4449
4450; plot data
4451   if keyword_set(timeseries) then  begin ; plot timeseries if keyword is set
4452         plottimeseries,stash
4453        endif else begin         ; otherwise plot points on a map
4454      plotpoints,stash
4455   endelse
4456
4457   prend2,view=0
4458
4459endif
4460
4461if (keyword_set(gif)) then begin
4462   thisDevice = !D.Name
4463
4464        Set_Plot, 'Z'         ; do graphics in the background
4465;  Device, Set_Resolution=[640,512], decomposed=0
4466;  Device, Set_Resolution=[800,512], decomposed=0
4467        Device, Set_Resolution=picsize, decomposed=0
4468        Erase                           ; clear any existing stuff
4469        !p.charsize=0.75
4470        setupct, r, g, b, coltable=stash.coltable, $
4471         white_on_black=stash.white_on_black    ; setup color table
4472
4473; plot data
4474   if keyword_set(timeseries) then  begin ; plot timeseries if keyword is set
4475         plottimeseries,stash
4476        endif else begin         ; otherwise plot points on a map
4477      plotpoints,stash
4478   endelse
4479
4480;  if (keyword_set(gif)) then begin
4481;        IMAGE_GRAB, /gif, filename=stash.outfile+'.gif',/overwrite       
4482   snapshot = TVRD()
4483        WRITE_GIF,stash.outfile+'.gif',snapshot, r, g, b
4484;        endif
4485        Device, Z_Buffer=1    ; reset graphics mode
4486        Set_Plot, thisDevice
4487        !p.charsize=0.0
4488
4489endif
4490
4491if (keyword_set(txt)) then begin
4492   if keyword_set(timeseries) then  begin ; plot timeseries if keyword is set
4493         plottimeseries,stash
4494        endif else begin         ; otherwise plot points on a map
4495      plotpoints,stash
4496;     print,'no txt version of plotpoints!'
4497   endelse
4498endif
4499
4500if (keyword_set(netcdf)) then begin
4501   if keyword_set(timeseries) then begin
4502         print,'no netcdf version of timeseries!'
4503        endif else begin
4504      plotpoints,stash       
4505        endelse
4506endif
4507
4508if (stash.filterout ne '') then begin
4509
4510   selpoints, stash, lonsel, latsel, innovsel, qcsel, daysel, obstypsel, obnumsel, numsel, typestr
4511
4512endif
4513
4514endelse            ; batch
4515 
4516END
4517
Note: See TracBrowser for help on using the repository browser.