;+
; NAME:
; thm_gmag_stackplot.pro
;
; PURPOSE:
; To create 3 PNG files displaying the H,D and Z components of the magnetic field
; from multiple GBO stations out of GMAG data that is stored in CDF file.
;
; CALLING SEQUENCE:
; thm_gmag_stackplot, date, duration, stack_shift=stack_shift, no_expose=no_expose, make_png=make_png
;
; INPUTS:
; date: The start of the time interval to be plotted. (Format: 'YYYY-MM-DD/hh:mm:ss')
; duration: The length of the interval being plotted. (Floating point number of days -> 12hr=0.5), default=1
; stack_shift: Space between stations on the y-axis (units are nanotesla), default=50
; no_expose: Set this keyword to prevent the plot from being printed to the screen.
; make_png: Set this keyword to make the 3 PNG files.
; max_deviation: Large spikes in the data (probably gliches) can screw up the y-axis scales. This keyword allows
; you to set the maximum deviation the data can go from the median; points that exceed this value are omitted.
; The default value is 500 nT.
; no_data_load: This keyword prevents new data from being loaded; the routine will try to plot existing data if it exists.
; hi_lat: Set this keyword to plot high latitude stations (above 49 degrees)
; lo_lat: Set this keyword to plot low latitude stations (below 49 below)
;
; OUTPUTS:
; None, but it creates 3 PNG files in the directory that IDL is being run.
;
; PROCEDURE:
; Read in data from CDF files; plot the data using tplot.pro routines;
; make PNG files with makegif.pro routine
;
; EXAMPLE:
; thm_gmag_stackplot, '2006-11-11',1,/make_png, max_deviation=500, stack_shift=200.
;
; MODIFICATION HISTORY:
; Written by: Matt Davis
; October 23, 2006 Initial version
;
;NOTE: This program is still in development. Features to be added:
; -generalizing the routine
;
;-
pro thm_gmag_stackplot, date, duration, stack_shift=stack_shift, no_expose=no_expose, make_png=make_png, max_deviation=max_deviation,no_data_load=no_data_load,hi_lat=hi_lat,lo_lat=lo_lat
;________________________________________________________________________________________________
;find all gmag CDF files on given date and put data into tplot variables
;________________________________________________________________________________________________
if not keyword_set(date) then begin
print, 'You must specify a date. (Format : YYYY-MM-DD/HH:MM:SS)'
return
endif
if not keyword_set(duration) then duration=1.
start_time=time_double(date)
end_time=start_time+86400.*duration
if not keyword_set(no_data_load) then thm_load_gmag, trange = [start_time,end_time],/verbose
;________________________________________________________________________________________________
;check that there is some valid data
;________________________________________________________________________________________________
tplotvars=tnames('thg_mag_????')
if tplotvars(0) eq '' then begin
print, 'Appropriate tplot variables could not be found!'
print, 'Searched for "thg_mag_????"'
print, 'Stackplot program was aborted.'
return
endif
;________________________________________________________________________________________________
;this section is in place to make sure that we are only looking at data in the time range specified,
;previously loaded gmag data from different time range is ignored
;________________________________________________________________________________________________
stations=['filler']
for i=0,n_elements(tplotvars)-1 do begin
get_data,tplotvars(i),data=dd
index_time=where(dd.x ge start_time and dd.x le end_time)
if index_time(0) ne -1 then stations=[stations,strmid(tplotvars(i),8,4)]
endfor
stations=stations(1:*)
;________________________________________________________________________________________________
;sort stations by longitude
;________________________________________________________________________________________________
lats=fltarr(n_elements(stations))
lons=fltarr(n_elements(stations))
for i=0,n_elements(stations)-1. do begin
get_data,'thg_mag_'+stations(i),dlimits=dl
lats(i)=float(dl.cdf.vatt.station_latitude)
lons(i)=float(dl.cdf.vatt.station_longitude)
endfor
hi_index=where(lats ge 49)
lo_index=where(lats lt 49)
hi_stat=stations(hi_index)
lo_stat=stations(lo_index)
hi_lon=lons(hi_index)
lo_lon=lons(lo_index)
hi_lat_lon_index=reverse(sort(hi_lon))
lo_lat_lon_index=reverse(sort(lo_lon))
hi_stat=hi_stat(hi_lat_lon_index)
lo_stat=lo_stat(lo_lat_lon_index)
;________________________________________________________________________________________________
;manipulate data into tplot variables for "relative" stack plots
;________________________________________________________________________________________________
if keyword_set(hi_lat) then a=0 else a=1 ; do high lat
if keyword_set(lo_lat) then b=1 else b=0 ; do low lat
if a gt b then begin ; no keywords then do both
a=0
b=1
endif
if not keyword_set(stack_shift) then set_default_shift=1 else set_default_shift=0
for w=a,b do begin
if w eq 0 then stations=hi_stat else stations=lo_stat
;create tplot variables of components
num_elements=n_elements(stations)
h_axis_range=dblarr(2)
d_axis_range=dblarr(2)
z_axis_range=dblarr(2)
if set_default_shift eq 1 and w eq 0 then stack_shift=500.
if set_default_shift eq 1 and w eq 1 then stack_shift=50.
if not keyword_set(max_deviation) then max_deviation=1000.
for i=0,num_elements-1 do begin
get_data,'thg_mag_'+stations(i),data=dd
index_time=where(dd.x ge start_time and dd.x le end_time)
themedian=strcompress(string(median(dd.y(index_time,0),/even),format='(f10.1)'),/remove_all)
hdata=dd.y(index_time,0)-median(dd.y(index_time,0),/even)
indexh=where(abs(hdata) le max_deviation)
hdata=hdata(indexh)+stack_shift*i
store_data,'thg_mag_'+stations(i)+'_h_rel',data={x:dd.x(index_time(indexh)),y:hdata},$
limits={labels:[strupcase(stations(i))+'='+string(themedian)]}
themedian=strcompress(string(median(dd.y(index_time,1),/even),format='(f10.1)'),/remove_all)
ddata=dd.y(index_time,1)-median(dd.y(index_time,1),/even)
indexd=where(abs(ddata) le max_deviation)
ddata=ddata(indexd)+stack_shift*i
store_data,'thg_mag_'+stations(i)+'_d_rel',data={x:dd.x(index_time(indexd)),y:ddata},$
limits={labels:[strupcase(stations(i))+'='+string(themedian)]}
themedian=strcompress(string(median(dd.y(index_time,2),/even),format='(f10.1)'),/remove_all)
zdata=dd.y(index_time,2)-median(dd.y(index_time,2),/even)
indexz=where(abs(zdata) le max_deviation)
zdata=zdata(indexz)+stack_shift*i
store_data,'thg_mag_'+stations(i)+'_z_rel',data={x:dd.x(index_time(indexz)),y:zdata},$
limits={labels:[strupcase(stations(i))+'='+string(themedian)]}
if max(hdata) gt h_axis_range(1) then h_axis_range(1)=max(hdata)
if max(ddata) gt d_axis_range(1) then d_axis_range(1)=max(ddata)
if max(zdata) gt z_axis_range(1) then z_axis_range(1)=max(zdata)
if min(hdata) lt h_axis_range(0) then h_axis_range(0)=min(hdata)
if min(ddata) lt d_axis_range(0) then d_axis_range(0)=min(ddata)
if min(zdata) lt z_axis_range(0) then z_axis_range(0)=min(zdata)
endfor
;________________________________________________________________________________________________
;create buffer tplot variables
;the purpose of these buffers is to 'trick' the plot routines into formatting the stackplots
;the way we want, it can be more reliable (as in doing what we want rather than what we tell it)
;than explicitly setting the plot format
;it also allows for the easy placement of the 'Median (nT)' label
;________________________________________________________________________________________________
bth=dd.y(index_time,0)
bbh=dd.y(index_time,0)
btd=dd.y(index_time,1)
bbd=dd.y(index_time,1)
btz=dd.y(index_time,2)
bbz=dd.y(index_time,2)
bth(*)=h_axis_range(1)+stack_shift
bbh(*)=h_axis_range(0)-stack_shift
btd(*)=d_axis_range(1)+stack_shift
bbd(*)=d_axis_range(0)-stack_shift
btz(*)=z_axis_range(1)+stack_shift
bbz(*)=z_axis_range(0)-stack_shift
store_data,'BUFFER_TOP_H',data={x:dd.x(index_time),y:bth}
store_data,'BUFFER_BOT_H',data={x:dd.x(index_time),y:bbh}
store_data,'BUFFER_TOP_D',data={x:dd.x(index_time),y:btd}
store_data,'BUFFER_BOT_D',data={x:dd.x(index_time),y:bbd}
store_data,'BUFFER_TOP_Z',data={x:dd.x(index_time),y:btz}
store_data,'BUFFER_BOT_Z',data={x:dd.x(index_time),y:bbz}
options,'BUFFER_TOP_H','ytitle',''
options,'BUFFER_BOT_H','ytitle',''
options,'BUFFER_TOP_D','ytitle',''
options,'BUFFER_BOT_D','ytitle',''
options,'BUFFER_TOP_Z','ytitle',''
options,'BUFFER_BOT_Z','ytitle',''
options,'BUFFER_TOP_H','labels',['Median (nT)']
options,'BUFFER_TOP_D','labels',['Median (nT)']
options,'BUFFER_TOP_Z','labels',['Median (nT)']
;________________________________________________________________________________________________
;Set various plotting options
;________________________________________________________________________________________________
tplotvars_h=['filler']
tplotvars_d=['filler']
tplotvars_z=['filler']
for i=0,n_elements(stations)-1 do begin
tplotvars_h=[tplotvars_h,tnames('thg_mag_'+stations(i)+'_h_rel')]
tplotvars_d=[tplotvars_d,tnames('thg_mag_'+stations(i)+'_d_rel')]
tplotvars_z=[tplotvars_z,tnames('thg_mag_'+stations(i)+'_z_rel')]
endfor
store_data,'BH',data=[tplotvars_h(1:*),'BUFFER_TOP_H','BUFFER_BOT_H']
store_data,'BD',data=[tplotvars_d(1:*),'BUFFER_TOP_D','BUFFER_BOT_D']
store_data,'BZ',data=[tplotvars_z(1:*),'BUFFER_TOP_Z','BUFFER_BOT_Z']
options,'BH','ytickformat','(a1)' ; this gets rid of the y-axis numbering
options,'BD','ytickformat','(a1)' ; this gets rid of the y-axis numbering
options,'BZ','ytickformat','(a1)' ; this gets rid of the y-axis numbering
num_hticks=(h_axis_range(1)+stack_shift)/100. - (h_axis_range(0)-stack_shift)/100.
h_factor=(byte(num_hticks/30.)+1)
num_hticks=byte(num_hticks/h_factor)-1
num_dticks=(d_axis_range(1)+stack_shift)/100. - (d_axis_range(0)-stack_shift)/100.
d_factor=(byte(num_dticks/30.)+1)
num_dticks=byte(num_dticks/d_factor)-1
num_zticks=(z_axis_range(1)+stack_shift)/100. - (z_axis_range(0)-stack_shift)/100.
z_factor=(byte(num_zticks/30.)+1)
num_zticks=byte(num_zticks/z_factor)-1
options,'BH','yticks',num_hticks
options,'BD','yticks',num_dticks
options,'BZ','yticks',num_zticks
hy_pos_ticks=findgen(byte( (h_axis_range(1)+stack_shift)/(100.*h_factor) +1.))*100.*h_factor
dy_pos_ticks=findgen(byte( (d_axis_range(1)+stack_shift)/(100.*d_factor) +1.))*100.*d_factor
zy_pos_ticks=findgen(byte( (z_axis_range(1)+stack_shift)/(100.*z_factor) +1.))*100.*z_factor
hy_neg_ticks=-100.*h_factor*reverse(findgen(byte( abs(h_axis_range(0)-stack_shift)/(100.*h_factor) +1.)))
dy_neg_ticks=-100.*d_factor*reverse(findgen(byte( abs(d_axis_range(0)-stack_shift)/(100.*d_factor) +1.)))
zy_neg_ticks=-100.*z_factor*reverse(findgen(byte( abs(z_axis_range(0)-stack_shift)/(100.*z_factor) +1.)))
if n_elements(hy_pos_ticks) gt 1 then hy_tickvalues=[hy_neg_ticks,hy_pos_ticks(1:*)] else hy_tickvalues=[hy_neg_ticks]
if n_elements(dy_pos_ticks) gt 1 then dy_tickvalues=[dy_neg_ticks,dy_pos_ticks(1:*)] else dy_tickvalues=[dy_neg_ticks]
if n_elements(zy_pos_ticks) gt 1 then zy_tickvalues=[zy_neg_ticks,zy_pos_ticks(1:*)] else zy_tickvalues=[zy_neg_ticks]
options,'BH','ytickv',hy_tickvalues
options,'BD','ytickv',dy_tickvalues
options,'BZ','ytickv',zy_tickvalues
options,'BH','yminor',10
options,'BD','yminor',10
options,'BZ','yminor',10
options,'BH','ytitle','Scale='+strcompress(string(100*h_factor,format='(f10.0)'))+'nT/major, '+strcompress(string(10*h_factor,format='(f10.0)'))+'nT/minor tickmark'
options,'BD','ytitle','Scale='+strcompress(string(100*d_factor,format='(f10.0)'))+'nT/major, '+strcompress(string(10*d_factor,format='(f10.0)'))+'nT/minor tickmark'
options,'BZ','ytitle','Scale='+strcompress(string(100*z_factor,format='(f10.0)'))+'nT/major, '+strcompress(string(10*z_factor,format='(f10.0)'))+'nT/minor tickmark'
;________________________________________________________________________________________________
;make plots/PNG files
;________________________________________________________________________________________________
If keyword_set(make_png) Then Begin
set_plot,'z'
device, set_resolution = [800, 800]
Endif Else Begin
set_plot,'x'
if keyword_set(no_expose) then window,0,xsize=800,ysize=800, /pixmap else window,0,xsize=800,ysize=800
Endelse
tplot_options,'region',[0.00,0.00,0.95,1.]
!p.color=0
!p.background=255
!p.charsize=1.2
time_stamp,/off
date_compress=strcompress(strmid(date,0,4)+strmid(date,5,2)+strmid(date,8,2),/remove_all)
timespan,date,duration
;--------------------------
if w eq 0 then begin
thetitle=string('Ground Magnetometer Data : High Latitude : H Component')
tplot,'BH',title=thetitle
if keyword_set(make_png) then makepng,'thg_l2_mag_hhla_'+string(date_compress)+'_v01',/no_expose
thetitle=string('Ground Magnetometer Data : High Latitude : D Component')
tplot,'BD',title=thetitle
if keyword_set(make_png) then makepng,'thg_l2_mag_dhla_'+string(date_compress)+'_v01',/no_expose
thetitle=string('Ground Magnetometer Data : High Latitude : Z Component')
tplot,'BZ',title=thetitle
if keyword_set(make_png) then makepng,'thg_l2_mag_zhla_'+string(date_compress)+'_v01',/no_expose
endif else begin
thetitle=string('Ground Magnetometer Data : Low Latitude : H Component')
tplot,'BH',title=thetitle
if keyword_set(make_png) then makepng,'thg_l2_mag_hlla_'+string(date_compress)+'_v01',/no_expose
thetitle=string('Ground Magnetometer Data : Low Latitude : D Component')
tplot,'BD',title=thetitle
if keyword_set(make_png) then makepng,'thg_l2_mag_dlla_'+string(date_compress)+'_v01',/no_expose
thetitle=string('Ground Magnetometer Data : Low Latitude : Z Component')
tplot,'BZ',title=thetitle
if keyword_set(make_png) then makepng,'thg_l2_mag_zlla_'+string(date_compress)+'_v01',/no_expose
endelse
;--------------------------
endfor ; w
end