;$Author: kenb-mac $
;$Date: 2007-01-24 14:23:38 -0800 (Wed, 24 Jan 2007) $
;$Header: /home/cdaweb/dev/control/RCS/plot_timeseries.pro,v 1.94 2006/12/11 16:47:41 kovalick Exp kovalick $
;$Locker: kovalick $
;$Revision: 225 $
;+------------------------------------------------------------------------
; NAME: PLOT_TIMESERIES
; PURPOSE: To generate a time series plot given the anonymous structures
; returned by the read_mycdf function.
; CALLING SEQUENCE:
; out = read_mycdf(Xstruct,Ystruct)
; INPUTS:
; Xstruct = structure containing the Epoch variable structure of the
; type returned by the read_mycdf structure.
; Ystruct = structure containing the variable to be plotted against
; the Epoch variable in the Xstruct parameter
; KEYWORD PARAMETERS:
; TSTART = Forces the time axis to begin at this Epoch time.
; TSTOP = Forces the time axis to end at this Epoch time.
; ELEMENTS = if set, then only these elements of a dimensional variable
; will be plotted.
; POSITION = If set, this routine will draw the plot(s) at this position
; of an existing window, rather than open a new one.
; FIRSTPLOT= Use this key in conjunction with the position keyword. Set
; this flag to indicate that the variable is the first in the
; window.
; LASTPLOT = Use this key in conjunction with the position keyword. Set
; this flag to indicate that the variable is the last in the
; window.
; PANEL_HEIGHT = vertical height, in pixels, of each panel
; CDAWEB = If set, the plot will have sufficient margin along the Z
; axis to hold a colorbar.
; GIF = If set, the plot will be a .gif file instead of Xwindow
; XSIZE = if set, forces the plot window to this width
; YSIZE = if set, forces the plot window to this height
; AUTO = if set, turns auto-scaling on
; NOGAPS = if set, eliminates data gap scanning
; NONOISE = if set, filter out values outside 3 sigma from mean
; IGNORE_DISPLAY_TYPE = if set, causes the attribute display_type to
; be ignored.
; NOSUBTITLE = if set, will not print 'time range = ' subtitle even after
; the last graph. Needed for timetext case. RCJ
; ONLYLABEL = if set, graph position is calculated but graph is not
; plotted. However, the x-axis label *is* plotted.
; Utilized by timetext case. RCJ
; SCATTER = if set, display a scatter plot (each point is plotted as a
; dot, no lines connect the dots. added on 5/14/2001 TJK
; COMBINE = if set, need to add the dataset name to the y axis label
; added 10/14/2003 - TJK.
; DEBUG = if set, turns on additional debug output.
; OUTPUTS:
; out = status flag, 0=0k, -1 = problem occured.
; AUTHOR:
; Richard Burley, NASA/GSFC/Code 632.0, Feb 22, 1996
; burley@nssdca.gsfc.nasa.gov (301)286-2864
; MODIFICATION HISTORY:
; 8/13/96 : R. Burley : Add NONOISE capability
;-------------------------------------------------------------------------
FUNCTION plot_timeseries, Xvar, Yvar, $
TSTART=TSTART,TSTOP=TSTOP,ELEMENTS=ELEMENTS,$
POSITION=POSITION,PANEL_HEIGHT=PANEL_HEIGHT,$
FIRSTPLOT=FIRSTPLOT,LASTPLOT=LASTPLOT,$
CDAWEB=CDAWEB,GIF=GIF,NOSUBTITLE=NOSUBTITLE,$
XSIZE=XSIZE,YSIZE=YSIZE, ONLYLABEL=ONLYLABEL,$
AUTO=AUTO,NOGAPS=NOGAPS,NOVALIDS,$
err_plus=err_plus,err_minus=err_minus, $
IGNORE_DISPLAY_TYPE=IGNORE_DISPLAY_TYPE,$
NONOISE=NONOISE,DEBUG=DEBUG,REPORT=REPORT,$
SCATTER=SCATTER,COMBINE=COMBINE,_EXTRA=EXTRAS
; Open report file if keyword is set
status = 0
;TJK added scatter plot capability on 5/14/2001
; RCJ 02/25/2005 Moved this to here because if 1 point then psym=4 (see
; somewhere below)
if keyword_set(SCATTER) then psym = 3 else psym = 0
if keyword_set(REPORT) then begin & reportflag=1L
a=size(REPORT) & if (a(n_elements(a)-2) eq 7) then $
OPENW,1,REPORT,132,WIDTH=132
endif else reportflag=0L
; Verify that both Xvar and Yvar are present
if (n_params() ne 2) then begin
print,'ERROR=Missing parameter to plot_timeseries function' & return,-1
endif
; Verify the type of the first parameter and retrieve the data
a = size(Xvar)
if (a(n_elements(a)-2) ne 8) then begin
print,'ERROR=1st parameter to plot_timeseries not a structure' & return,-1
endif else begin
a = tagindex('DAT',tag_names(Xvar))
if (a(0) ne -1) then times = Xvar.DAT $
else begin
a = tagindex('HANDLE',tag_names(Xvar))
if (a(0) ne -1) then handle_value,Xvar.HANDLE,times $
else begin
print,'ERROR=1st parameter does not have DAT or HANDLE tag' & return,-1
endelse
b = size(times)
; help, /struct, b
;TJK 6/27/2006 - added check for "9" - cdf_epoch16, "5" is cdf_epoch
if ((b(n_elements(b)-2) eq 5) or (b(n_elements(b)-2) eq 9)) then begin
;print, 'Epoch parameter ok'
endif else begin
print,'ERROR=1st parameter datatype not CDF EPOCH' & return,-1
endelse
endelse
endelse
tszck=size(times)
if(tszck(tszck(0)+2) ne 1) then $ ; RTB added to prevent reform(scalar)
times = reform(times) ; eliminate any redundant dimensions
; Verify the type of the second parameter and retrieve the data
a = size(Yvar)
if (a(n_elements(a)-2) ne 8) then begin
print,'ERROR=2nd parameter to plot_timeseries not a structure' & return,-1
endif else begin
YTAGS = tag_names(Yvar) ; avoid multiple calls to tag_names
a = tagindex('DAT',YTAGS)
if (a(0) ne -1) then THEDATA = Yvar.DAT $
else begin
a = tagindex('HANDLE',YTAGS)
if (a(0) ne -1) then handle_value,Yvar.HANDLE,THEDATA $
else begin
print,'ERROR=2nd parameter does not have DAT or HANDLE tag' & return,-1
endelse
endelse
endelse
szck=size(thedata)
if(szck(szck(0)+2) ne 1) then $ ; RTB added to prevent reform(scalar)
thedata = reform(thedata) ; eliminate any redundant dimensions
; Verify type of data and determine the number of panels that will be plotted
; and which elements of the data array are to be plotted.
a = size(thedata) & b = a(n_elements(a)-2) & thedata_size = a
;TJK changed the following line on 5/26/2000 so that unisigned datatypes are
;also acceptable.
;if ((b eq 0) OR (b gt 5)) then begin
if ((b eq 0) OR (b gt 6 and b lt 12)) then begin
print,'STATUS=datatype indicates that data is not plottable' & return,-1
endif else begin
case a(0) of
0 : begin
; what follows is a little trick to allow us to plot
; one single point. RCJ 02/2001
thedata=[thedata,thedata]
a = size(thedata)
b = a(n_elements(a)-2) & thedata_size = a
times=[times,times]
psym=4
symsize=2
;print,'ERROR=Single data points are not plottable' & return,-1
;print,'STATUS=Re-select longer time interval. Single data points are not plottable'
print,'STATUS: Found one single point.....'
num_panels = 1L & elist=0
;return,-1
end
1 : begin
num_panels = 1L & elist=0
;TJK 2/10/2004 - put this check in, so that we're not trying to plot unequal arrays
;print, n_elements(thedata), n_elements(times)
if (n_elements(thedata) ne n_elements(times)) then begin
; print,'ERROR=Single data points are not plottable'
print,'STATUS=Re-select longer time interval. Arrays of values for a single time are not plottable'
return, -1
endif
end
2 : begin ; #panels determined by dimensionality or by display type
elist=indgen(a(1))
;TJK change to check if elements is defined
; if keyword_set(ELEMENTS) then elist = ELEMENTS $
if (n_elements(ELEMENTS) gt 0) then elist = ELEMENTS $
else begin
if NOT keyword_set(IGNORE_DISPLAY_TYPE) then begin
b = tagindex('DISPLAY_TYPE',YTAGS)
if (b(0) ne -1) then begin ; evaluate the display type
c = strupcase(Yvar.(b(0))) & c = break_mystring(c,delimiter='>')
if ((c(0) eq 'TIME_SERIES')AND(n_elements(c) gt 1)) then begin
d = break_mystring(c(1),delimiter=',')
elist = long(d) & elist = elist -1
endif
endif
endif
endelse
end
else: begin
print,'ERROR=Cannot plot data with > 2 dimensions' & return,-1
end
endcase
endelse
num_panels = n_elements(elist)
;print, 'num_panels set to n_elements(elist) = ',num_panels
; Determine the proper start and stop times of the plot
tbegin = times(0) & tend = times(n_elements(times)-1) ; default to data
if keyword_set(TSTART) then begin ; set tbegin
tbegin = TSTART & tbegin16 = TSTART & a = size(TSTART)
if (a(n_elements(a)-2) eq 7) then begin ;if tstart is a string, convert it
tbegin = encode_CDFEPOCH(TSTART)
tbegin16 = encode_CDFEPOCH(TSTART,/EPOCH16);TJK added for use when data is epoch16
endif
endif
if keyword_set(TSTOP) then begin ; set tend
tend = TSTOP & tend16 = TSTOP & a = size(TSTOP)
if (a(n_elements(a)-2) eq 7) then begin ;if tstop is a string, convert it
tend = encode_CDFEPOCH(TSTOP)
tend16 = encode_CDFEPOCH(TSTOP,/EPOCH16);TJK added for use when data is epoch16
endif
endif
; Compare the range of the time data to the requested start and stop times
pad_front = 0L & pad_end = 0L
;print, tbegin, tend
;print, '1st time val from data ',times(0)
;print, 'last time val from data ',times(n_elements(times)-1)
;TJK 7/20/2006 test to see if time values are standard epoch or
;epoch16, set flag and change the tend and tbegin values to their
;epoch16 counterparts so we don't have to change so much code down
;below
ep16 = 0
if (size(times(0),/tname) eq 'DCOMPLEX')then begin
ep16 = 1
tend = tend16
tbegin = tbegin16
;print, 'data is dcomplex, ',tend, ' and ',tbegin
endif
if (!version.release ge '6.2' and ep16) then begin
if (cdf_epoch_compare(times(0), tbegin)) then begin
if keyword_set(DEBUG) then print,'Padding front of times...'
times = [tbegin,times] & pad_front = 1L
endif
if (cdf_epoch_compare(tend, times(n_elements(times)-1))) then begin
if keyword_set(DEBUG) then print,'Padding end of times...'
times = [times,tend] & pad_end = 1L
endif
endif else begin
if (tbegin lt times(0)) then begin
if keyword_set(DEBUG) then print,'Padding front of times...'
times = [tbegin,times] & pad_front = 1L
endif
if (tend gt times(n_elements(times)-1)) then begin
if keyword_set(DEBUG) then print,'Padding end of times...'
times = [times,tend] & pad_end = 1L
endif
endelse
; Determine the first and last data time points to be plotted
rbegin = 0L & w = where(times ge tbegin,wc)
if (wc gt 0) then rbegin = w(0)
;print,rbegin
rend = n_elements(times)-1 & w = where(times le tend,wc)
if (wc gt 0) then rend = w(n_elements(w)-1)
if (rbegin ge rend) then begin
print,'STATUS=No data within specified time range.' ;& return,-1
endif
if not (keyword_set(nosubtitle)) then begin
; Create a subtitle for the plots showing the data start and stop times
CDF_EPOCH,tbegin,year,month,day,hour,minute,second,milli,/BREAK
ical,year,doy,month,day,/idoy
subtitle = 'TIME RANGE='+strtrim(string(year),2)+'/'+strtrim(string(month),2)
subtitle = subtitle + '/' + strtrim(string(day),2)
subtitle = subtitle + ' (' + strtrim(string(doy),2) + ') to '
CDF_EPOCH,tend,year,month,day,hour,minute,second,milli,/BREAK
ical,year,doy,month,day,/idoy
subtitle = subtitle + strtrim(string(year),2)+'/'+strtrim(string(month),2)
subtitle = subtitle + '/' + strtrim(string(day),2)
subtitle = subtitle + ' (' + strtrim(string(doy),2) + ')'
endif else subtitle=''
; Convert the time array into seconds since tbegin
CDF_EPOCH,tbegin,year,month,day,hour,minute,second,milli,/BREAK
CDF_EPOCH,a,year,month,day,0,0,0,0,/COMPUTE_EPOCH
temp = make_array(n_elements(times),/dcomplex)
;TJK 7/20/2006 - new section to "do the math" on the epoch16 values
if (ep16)then begin
for i = 0L, n_elements(times)-1 do begin
; double in seconds from first of day
ep_diff = cdf_epoch_diff (times(i), a)
; help, ep_diff
; print, 'diff = ', ep_diff
temp(i) = ep_diff/1000.d0
endfor
times = temp
endif else begin
times = (times - a) / 1000.d0 ; double in seconds from first of day
endelse
;print,'after conversion ', times
julday = ymd2jd(year,month,day)
; Determine label for time axis based on time range
xranger = lonarr(2)
;TJK 7/18/2006 - if ep16, replace w/ the "real parts" of the 1st and last element
;of the times array, which was already "converted" (reduced by a and /1000) above.
if (ep16) then begin
xranger(0) = real_part(times(0))
xranger(1) = real_part(times(n_elements(times)-1))
endif else begin
xranger(0) = (tbegin-a)/1000
xranger(1) = (tend-a)/1000
endelse
trange = xranger(1) - xranger(0)
if (trange le 60.0) then tform='h$:m$:s$.f$@y$ n$ d$' $
else tform='h$:m$:s$@y$ n$ d$'
; Determine if a new window should be opened for the plots
new_window = 1L ; initialize assuming a new window needs to be created
if keyword_set(POSITION) then begin ; adding to existing plot
a = size(POSITION) & b = n_elements(a)
if ((a(b-1) ne 4)OR(a(b-2) ne 3)) then begin
print,'ERROR=Invalid value for POSITION keyword' & return,-1
endif
if keyword_set(PANEL_HEIGHT) then begin ; verify it
a = size(PANEL_HEIGHT) & b = n_elements(a)
if ((a(b-2) le 1)OR(a(b-2) gt 5)) then begin
print,'ERROR=Invalid value for PANEL_HEIGHT keyword' & return,-1
endif else psize = PANEL_HEIGHT
endif else begin
print,'ERROR=PANEL_HEIGHT keyword must be specified with POSITION keyword'
return,-1
endelse
if keyword_set(FIRSTPLOT) then clear_plot = 0L else clear_plot = 1L
new_window = 0L ; no new window needed
endif
; Determine the size of the new window to be created
if (new_window eq 1) then begin
if keyword_set(GIF) then begin
xs = 640 & ys = 512 & psize = 100 ; set default gif sizes
if keyword_set(XSIZE) then xs = XSIZE ; override if keyword present
if keyword_set(YSIZE) then ys = YSIZE ; override if keyword present
if keyword_set(PANEL_HEIGHT) then begin
psize = PANEL_HEIGHT & ys = (psize * num_panels) + 100
endif else psize = ((ys-100) / num_panels)
endif else begin ; generating an X-window
a = lonarr(2) & DEVICE,GET_SCREEN_SIZE=a ; get device resolution
xs = (a(0)*0.66) & ys = (a(1)*0.66) ; compute defaults
if keyword_set(XSIZE) then xs = XSIZE ; override if keyword present
if keyword_set(YSIZE) then ys = YSIZE ; override if keyword present
if keyword_set(PANEL_HEIGHT) then begin
psize = PANEL_HEIGHT
ys = (psize * num_panels) + 100
if (ys gt a(1)) then begin
print,'ERROR=Computed window Ysize greater than device resolution'
return,-1
endif
endif else psize = ((ys-100) / num_panels)
endelse
if (psize lt 50) then begin ; sanity check for #pixels per panel
print,'ERROR=Insufficient resolution for a ',num_panels,' panel plot'
return,-1
endif
endif
; Initialize plotting position arrays and flags
if keyword_set(POSITION) then ppos = POSITION $
else begin
ppos = fltarr(4) ; create position array
ppos(0) = 100 ; default plot x origin
ppos(2) = (xs - 40) ; default plot x corner
ppos(1) = (ys - 30) - psize ; 1st plot y origin
ppos(3) = (ys - 30) ; 1st plot y corner
if keyword_set(CDAWEB) then ppos(2) = xs - 100 ; set margin for spectrogram
endelse
; Determine the title for the window or gif file
if (new_window eq 1) then begin
a = tagindex('SOURCE_NAME',YTAGS)
if (a(0) ne -1) then b = Yvar.SOURCE_NAME else b = ''
a = tagindex('DESCRIPTOR',YTAGS)
if (a(0) ne -1) then b = b + ' ' + Yvar.DESCRIPTOR
window_title = b
endif
; Output debug output if requested
;if keyword_set(DEBUG) then begin
; print,'TIMEAXIS FORMAT=',tform
; print,'PANEL HEIGHT=',psize
; print,'#of PANELS=',num_panels
; print,'ELEMENT LIST=',elist+1
; print,'pad_front=',pad_front,' pad_end=',pad_end
;endif
; Create the new window or the gif file
if (new_window eq 1) then begin
if keyword_set(GIF) then begin
a = size(GIF) & if (a(n_elements(a)-2) ne 7) then GIF = 'idl.gif'
deviceopen,6,fileOutput=GIF,sizeWindow=[xs,ys]
endif else begin ; open x-window display
window,/FREE,XSIZE=xs,YSIZE=ys,TITLE=window_title
clear_plot = 0L ; initialize clear plot flag
endelse
endif
; Determine the fill value for the Y data and valid min and valid max values
a = tagindex('FILLVAL',YTAGS)
Yfillval = 1.0e31
if (a(0) ne -1) then begin
if (Yvar.FILLVAL ne '') then Yfillval = Yvar.FILLVAL
endif
if keyword_set(err_plus) and keyword_set(err_minus) then begin
if (n_elements(err_plus) ne n_elements(thedata) or $
n_elements(err_minus) ne n_elements(thedata)) then begin
err_plus=0
err_minus=0
print,'Plot_timeseries: Could not plot error bars'
endif
endif
; EXTRACT THE DATA FOR EACH PANEL AND PLOT
for i=0,num_panels-1 do begin
; extract data for a single panel from the data array
if (thedata_size(0) eq 1) then begin
mydata = thedata
if keyword_set (err_plus) then myerr_plus = err_plus
if keyword_set (err_minus) then myerr_minus = err_minus
endif else begin
mydata = thedata((elist(i)),*)
if keyword_set (err_plus) then myerr_plus = err_plus((elist(i)),*)
if keyword_set (err_minus) then myerr_minus = err_minus((elist(i)),*)
endelse
mydata = reform(mydata) ; remove any extraneous dimensions
if keyword_set (err_plus) then myerr_plus=reform(myerr_plus)
if keyword_set (err_minus) then myerr_minus=reform(myerr_minus)
; pad the beginning and end of data if extra time points were added
if (pad_front) then begin
mydata = [Yfillval,mydata] ; add fill point to front
if keyword_set (err_plus) then myerr_plus=[myerr_plus(0),myerr_plus]
if keyword_set (err_minus) then myerr_minus=[myerr_minus(0),myerr_minus]
endif
if (pad_end) then begin
mydata = [mydata,Yfillval] ; add fill point to back
if keyword_set (err_plus) then myerr_plus=[myerr_plus,myerr_plus(n_elements(myerr_plus)-1)]
if keyword_set (err_minus) then myerr_minus=[myerr_minus,myerr_minus(n_elements(myerr_minus)-1)]
endif
; screen out data points which are outside the plotting time range
; Check data before plotting
rrend=n_elements(mydata)
; if(rrend lt rend) then begin
; print, "STATUS=No Data Available"
; return, -1
; endif
;TJK modified this so that it will just produce a blank plot
;in the case where there's data for some variables in this
;dataset but not all. It was too hard to figure out how to
;make the s/w "skip" a variable since the size of the window
;and whatnot is done in plotmaster w/ no knowledge of the data found
;for each variable.
;TJK added 1/29/99 setting of flag
nogood = 0; FALSE
if(rrend lt rend) then rend = rrend-1
mydata = mydata(rbegin:rend)
if keyword_set (err_plus) then myerr_plus = myerr_plus(rbegin:rend)
if keyword_set (err_minus) then myerr_minus = myerr_minus(rbegin:rend)
mytimes = times(rbegin:rend)
; screen out fill data from the data and the times array
w = where(mydata ne Yfillval,non_fillcount)
;TJK 8/5/2004 - add this check for inifinite or NaN "values" since
;we now have datasets that have whole arrays of NaN, e.g. wi_phsp_3dp
n_goodvals = 0 ; need to initialize
if (non_fillcount ne 0) then n = where(finite(mydata(w)) eq 1,n_goodvals)
;
if (non_fillcount ne 0 and n_goodvals gt 0) then begin
mydata = mydata(w) & mytimes = mytimes(w)
if keyword_set (err_plus) then myerr_plus=myerr_plus(w)
if keyword_set (err_minus) then myerr_minus=myerr_minus(w)
w=0
endif else begin
w=0
nogood = 1
endelse
; screen out data outside validmin and validmax values
if ((NOT keyword_set(NOVALIDS))AND(non_fillcount gt 0)and $
(n_goodvals gt 0)) then begin
; determine validmin and validmax values
a = tagindex('VALIDMIN',YTAGS)
if (a(0) ne -1) then begin & b=size(Yvar.VALIDMIN)
if (b(0) eq 0) then Yvmin = Yvar.VALIDMIN $
else Yvmin = Yvar.VALIDMIN(elist(i))
endif else Yvmin = 1.0e31
a = tagindex('VALIDMAX',YTAGS)
if (a(0) ne -1) then begin & b=size(Yvar.VALIDMAX)
if (b(0) eq 0) then Yvmax = Yvar.VALIDMAX $
else Yvmax = Yvar.VALIDMAX(elist(i))
endif else Yvmax = 1.0e31
; proceed with screening
w = where(((mydata gt Yvmax)OR(mydata lt Yvmin)),wc)
if (wc gt 0) then begin
if keyword_set(DEBUG) then print,wc,' values outside VALIDMIN/MAX'
w = where(((mydata le Yvmax)AND(mydata ge Yvmin)),wb)
if (wb gt 0) then begin
mydata=mydata(w) & mytimes=mytimes(w)
if keyword_set (err_plus) then myerr_plus=myerr_plus(w)
if keyword_set (err_minus) then myerr_minus=myerr_minus(w)
endif else begin
a = tagindex('FIELDNAM',YTAGS)
if (a(0) ne -1) then ylabel = Yvar.(a(0))
print,'STATUS=No data for at least one component of ',ylabel,' variable.'
;TJK want to continue since this might be one of several panels - instead
;lets set a flag and check it below.
; return,0
nogood = 1; TRUE
endelse
endif
endif
; RCJ 10/16/02 If this is the first plot and the data is not good, at least create
; and invisible plot so subsequent graphs won't fall on a black graph area.
; In the case where this is the only graph requested
; at least the user won't end up looking at a black graph area.
; TJK 8/5/2004 added check for i eq 0 which means 1st of multiple panels
if (nogood and firstplot and (i eq 0)) then begin
plot,[0,1],[0,1],/nodata,ystyle=8+4,xstyle=8+4
endif
if (nogood eq 0) then begin
; screen out data outsize 3 standard deviations from the mean
if keyword_set(NONOISE) then begin
;semiMinMax,mydata,Sigmin,Sigmax
; RCJ 05/01/2006 Replaced call to semiminmax w/ call to three_sigma
sigminmax=three_sigma(mydata)
sigmin=sigminmax.(0)
sigmax=sigminmax.(1)
w = where(((mydata gt Sigmax)OR(mydata lt Sigmin)),wc)
if (wc gt 0) then begin
if keyword_set(DEBUG) then print,wc,' values outside 3-sigma...'
w = where(((mydata le Sigmax)AND(mydata ge Sigmin)),wb)
if (wb gt 0) then begin
mydata=mydata(w) & mytimes=mytimes(w)
if keyword_set (err_plus) then myerr_plus=myerr_plus(w)
if keyword_set (err_minus) then myerr_minus=myerr_minus(w)
endif
endif
endif
; determine the yaxis scale type, natural or logarithmic
yscaletype = 0L ; initialize assuming natural
a = tagindex('SCALETYP',YTAGS)
if (a(0) ne -1) then begin
if (strupcase(Yvar.SCALETYP) eq 'LOG') then yscaletype = 1L
endif
; screen non-positive data values if creating a logarithmic plot
if (yscaletype eq 1) then begin
wle = where(mydata le 0.0,wcle)
if (wcle gt 0) then begin
w = where(mydata gt 0.0,wc)
if (wc gt 0) then begin ;if there are good values
;TJK 10/01/2004 - change from just throwing out values <=0 to reassigning them
wmin = min(mydata(w));get smallest real value above zero
wmin = wmin/2 ; make it less than the real smallest value - TJK 10/22/2004
mydata(wle) = wmin ;assign all valid values zero and below to a really small value
;above zero
;SAVE wle and wcle for use in the plot command below
;TJK 10/01/2004 - remove this and do above mydata = mydata(w) & mytimes = mytimes(w)
w = where(mydata gt 0.0,wc)
if (wc gt 0) then begin ;if there are good values
if keyword_set (err_plus) then myerr_plus=myerr_plus(w)
if keyword_set (err_minus) then myerr_minus=myerr_minus(w)
endif
w=0
endif
endif
endif
; determine the proper scale for the Y axis
;TJK initialize to scales based on the real data. Because we're now enabling
; a "noauto" keyword in the time_series synax for display_type, we want to
; initialize ymin and max to more realistic values just in case the scalemin
; and max attributes aren't set.
;
goodvals = where(finite(mydata) eq 1, ngoodvals)
ymin = min(mydata(goodvals),MAX=ymax)
;TJK replaced w/ above ymax = 1.0 & ymin = 0.0 ; initialize to any value
a = tagindex('SCALEMIN',YTAGS)
if (a(0) ne -1) then begin & b=size(Yvar.SCALEMIN)
if (b(0) eq 0) then ymin = Yvar.SCALEMIN $
else ymin = Yvar.SCALEMIN(elist(i))
endif
a = tagindex('SCALEMAX',YTAGS)
if (a(0) ne -1) then begin & b=size(Yvar.SCALEMAX)
if (b(0) eq 0) then ymax = Yvar.SCALEMAX $
else ymax = Yvar.SCALEMAX(elist(i))
endif
if (keyword_set(AUTO)) then begin ; autoscale based on valid data values
if (non_fillcount gt 0) then begin ; cant autoscale if all fill data
ymax = 0.0
; replace this with next 2 lines to ignore NaN values ymin = min(mydata,MAX=ymax)
goodvals = where(finite(mydata) eq 1, ngoodvals)
ymin = min(mydata(goodvals),MAX=ymax)
endif
endif
; print, 'TJK debug - data ymin and max ',ymin, ymax
; quality check the y scales
if (ymax eq ymin) then begin
;ymax = ymax + 1
;ymin = ymin - 1
; RCJ 02/25/2005 10% of the value is probably better.
ymax = ymax + (ymax * .1)
ymin = ymin - (ymin * .1)
; RCJ 05/02/2006 If ymin=ymax=0. (case found after removing noise
; from AC_H1_EPM data) then we don't really want a log plot:
if ymax eq 0. and ymin eq 0. then begin
yscaletype=0
ymin=-.1 & ymax=.1
endif
endif
;print, 'DEBUG, ymin and max = ',ymin,ymax ;TJK DEBUG
yranger=[ymin,ymax]
; RCJ 11/15/02 Adjusting graph yrange according to error bars.
; The graphs were not improved if the y-axis was in log scale so no correction
; is made for those cases.
if (yscaletype ne 1) then begin
if keyword_set(err_minus) or keyword_set(err_plus) then $
print,'Adjusting ymin and ymax according to error bars...'
if keyword_set (err_minus) then begin
q=where(mydata eq ymin)
; if q ne -1, than ymin exists in mydata and I only need its index number.
; if q=-1, when ymin=ymax, then just take ymin
if q[0] ne -1 then yranger[0]= ymin-max(myerr_minus[q]) $
else yranger[0]= ymin
endif
if keyword_set (err_plus) then begin
q=where(mydata eq ymax)
if q[0] ne -1 then yranger[1] = ymax+max(myerr_plus[q]) $
else yranger[1] = ymax
endif
endif
if ((yscaletype eq 1)AND(yranger(0) le 0)) then yranger(0) = 0.00001
;print, 'DEBUG, yranger = ', yranger ;TJK DEBUG
; Determine the proper labeling for the y axis
ylabel = '' & yunits = '' & yds = '' ; initialize
if keyword_set(COMBINE) then begin
a = tagindex('LOGICAL_SOURCE',YTAGS)
if (a(0) ne -1) then yds = strupcase(Yvar.(a(0)))
endif
a = tagindex('FIELDNAM',YTAGS)
if (a(0) ne -1) then ylabel = Yvar.(a(0))
a = tagindex('LABLAXIS',YTAGS)
if (a(0) ne -1) then ylabel = Yvar.(a(0))
a = tagindex('LABL_PTR_1',YTAGS)
if (a(0) ne -1) then begin
if (Yvar.(a(0))(0) ne '') then ylabel = Yvar.(a(0))(elist(i))
endif
a = tagindex('UNITS',YTAGS)
if (a(0) ne -1) then yunits = Yvar.(a(0))
a = tagindex('UNIT_PTR',YTAGS)
if (a(0) ne -1) then begin
if (Yvar.(a(0))(0) ne '') then yunits = Yvar.(a(0))(elist(i))
endif
;
; should add a condition here: only do the next 2 lines if sscweb??????
;
coord=str_sep(yvar.varname,'_')
if (n_elements(coord) gt 1) and (coord(0) eq 'XYZ') then ylabel=ylabel+' ('+coord(1)+')'
if (n_elements(yds) gt 0) then begin
ylabel = yds + '!C' + ylabel + '!C' + yunits ;TJK removed for space 10/12/2006 + '!C'
endif else ylabel = ylabel + '!C' + yunits ;same as above + '!C'
; compare the size of the yaxis label to the panel height
ycsize = 1.0 & ylength = strlen(ylabel)
; if ((!d.y_ch_size * ylength) gt psize) then begin
; ratio = float(!d.y_ch_size * ylength) / float(psize)
; ycsize = 1.0 - (ratio/8.0) + 0.1
; endif
;TJK since this label is written vertically, you want to use the
;characters "x" size in determining the overall size that will fit.
if ((!d.x_ch_size * ylength) gt psize) then begin
ratio = float(!d.x_ch_size * ylength) / float(psize)
ycsize = 1.0 - (ratio/8.0) + 0.1
endif
; search for data gaps
if keyword_set(NOGAPS) then datagaps = -1 else datagaps = find_gaps(mytimes)
; produce debug output
if keyword_set(DEBUG) then begin
print,' Yscales=',yranger & print,' Yscaletype=',yscaletype
endif
; generate an empty plot frame
; rtb added 12/98
;!p.charsize=1.0
if keyword_set(onlylabel) then begin
plot,mytimes,mydata,XSTYLE=4+1,ystyle=4+1,/NODATA,$
XRANGE=xranger,POSITION=ppos,/DEVICE,NOERASE=clear_plot
!y.crange(0)=!y.crange(1)
timeaxis_text,JD=julday,form=tform,/onlylabel
goto, skipped_graph
endif else begin
;TJK 10/22/2004 - changed ystyle from 2 to 1 - 1 forces exact range which
;is what I thought we were doing long ago!
; plot,mytimes,mydata,/NODATA,YTITLE=ylabel,YRANGE=yranger,YSTYLE=2,$
;TJK 10/26/2004 - try something different to get better labeling
;
if (yscaletype) then begin ;log scaling -
;10/28/2004 - TJK - get the nice y range that IDL
;would use and set
;any values le 0.0 to the bottom scale value and
;call plot again forcing the exact scale ranges to
;the min and max returned in ytick_get
plot,mytimes,mydata,/NODATA,YTITLE=ylabel,YRANGE=yranger,YSTYLE=2+4,$
YLOG=yscaletype,XSTYLE=4+1,XRANGE=xranger,POSITION=ppos,/DEVICE,$
NOERASE=clear_plot,CHARSIZE=ycsize,_EXTRA=EXTRAS,ytick_get=yticks
;reassign the values originally le 0.0 to the scale min (returned by IDL)
if (wcle gt 0) then begin
mydata(wle) = yticks(0)
if keyword_set(DEBUG) then print,'Log scaling - reassigning values le 0 to ',yticks(0)
endif
if (n_elements(yticks) gt 0) then begin
yranger(0) = min(yticks)
yranger(1) = max(yticks)
endif
;use exact yscale values for log plots (ystyle=1).
plot,mytimes,mydata,/NODATA,YTITLE=ylabel,YRANGE=yranger,YSTYLE=1,$
YLOG=yscaletype,XSTYLE=4+1,XRANGE=xranger,POSITION=ppos,/DEVICE,$
NOERASE=clear_plot,CHARSIZE=ycsize,_EXTRA=EXTRAS
endif else begin ;linear scaling - use the IDL defined yscales (ystyle=2)
;TJK 7/17/2006 - check for values that IDL likes to label in exponential
;form - force them to just F notation, otherwise the variable labels
;get forced off the left side. Leave the really large numbers to IDL for now.
if (abs(yranger(1)) gt 999 and abs(yranger(1)) lt 99999) then begin
plot,mytimes,mydata,/NODATA,YTITLE=ylabel,YRANGE=yranger,YSTYLE=2,$
YLOG=yscaletype,XSTYLE=4+1,ytickformat='(F11.1)',XRANGE=xranger,POSITION=ppos,/DEVICE,$
NOERASE=clear_plot,CHARSIZE=ycsize,_EXTRA=EXTRAS
;TJK 7/28/2006 - check for small values where psize is ge 200 (our
;large panels) use a ytickformat that won't push the labels off the left side
endif else if ((abs(yranger(1)) lt 1.0) and (psize ge 200)) then begin
plot,mytimes,mydata,/NODATA,YTITLE=ylabel,YRANGE=yranger,YSTYLE=2,$
YLOG=yscaletype,XSTYLE=4+1,ytickformat='(F9.6)',XRANGE=xranger,POSITION=ppos,/DEVICE,$
NOERASE=clear_plot,CHARSIZE=ycsize,_EXTRA=EXTRAS
endif else begin
plot,mytimes,mydata,/NODATA,YTITLE=ylabel,YRANGE=yranger,YSTYLE=2,$
YLOG=yscaletype,XSTYLE=4+1,XRANGE=xranger,POSITION=ppos,/DEVICE,$
NOERASE=clear_plot,CHARSIZE=ycsize,_EXTRA=EXTRAS
endelse
endelse
;help,mydata,err_minus,err_plus
if keyword_set(err_plus) and keyword_set(err_minus) then $
cdaweb_errplot,mytimes,mydata-myerr_minus,mydata+myerr_plus
timeaxis_text,JD=julday,/NOLABELS,TICKLEN=-2.0
endelse
; if any plottable data exists then overplot the data into the frame
;TJK added scatter plot capability on 5/14/2001
;if keyword_set(SCATTER) then psym = 3 else psym = 0
if (non_fillcount ne 0) then begin
if (datagaps(0) eq -1) then oplot,mytimes,mydata,psym=psym,symsize=symsize $
else begin
start = 0L ; overplot each data segment
for j=0L,n_elements(datagaps)-1 do begin
stop = datagaps(j) ; get last element of segment
;oplot,mytimes(start:stop),mydata(start:stop)
oplot,mytimes(start:stop),mydata(start:stop), psym=psym
start = stop + 1; reset start element for next oplot
endfor
;oplot,mytimes(start:*),mydata(start:*) ; oplot last segment
oplot,mytimes(start:*),mydata(start:*), psym=psym ; oplot last segment
endelse
endif
;endif else print, "STATUS=No Data Available" & status = -2
; Adjust plot position and flags for next plot
ppos(3) = ppos(1) & ppos(1) = ppos(1) - psize & clear_plot=1
endif ;nogood
endfor
; RCJ 10/16/02 If this data is good and it's not the first plot then
; we are guaranteed not to have a black graph area and we can go ahead with
; the timeaxis_text calls.
; RCJ 11/18/02 It looks like 'and firstplot' is not needed. Will test
; in dev and bring it back if needed.
;if not (nogood and firstplot) then begin
if not (nogood) then begin
; draw the time axis by default or if lastplot flag is specified
if keyword_set(POSITION) then begin
if keyword_set(LASTPLOT) then begin
timeaxis_text,FORM=tform,JD=julday,title=subtitle,CHARSIZE=0.9
endif
endif else begin
timeaxis_text,FORM=tform,JD=julday,title=subtitle,CHARSIZE=0.9
if keyword_set(GIF) then begin
print,'not yet titling gifs from within plot_timeseries'
deviceclose
endif
endelse
endif
skipped_graph:
return,status
end