; ROC and STONE curve calculations for the SWMF-Dst comparison. close,1 device,decompose=0 monthmin=fltarr(11,12) monthmax=fltarr(11,12) monthdays=fltarr(11,12) first_year=2015. ; first day of values for each month: start at year 2007 monthmin(0,*)=[1,13,1,10,1,1,2,9,1,1,1,1] ; 2007 monthmin(1,*)=[1,6,1,11,3,16,1,4,8,1,1,1] ; 2008 monthmin(2,*)=[1,1,1,1,1,1,1,1,1,1,1,1] ; 2009 monthmin(3,*)=[1,1,1,1,1,1,1,1,1,1,1,1] ; 2010 monthmin(4,*)=[1,1,1,1,5,1,1,1,9,1,1,1] ; 2011 monthmin(5,*)=[9,1,1,1,1,1,2,1,1,1,1,1] ; 2012 monthmin(6,*)=[1,1,1,1,9,1,1,1,1,1,1,1] ; 2013 monthmin(7,*)=[1,1,1,1,1,1,1,1,1,1,4,1] ; 2014 monthmin(8,*)=[1,1,1,1,1,1,1,1,1,1,1,1] ; 2015 monthmin(9,*)=[1,1,1,1,1,1,1,1,1,1,1,1] ; 2016 monthmin(10,*)=[1,1,1,1,1,1,1,1,1,1,1,1] ; 2017 ; last day of values for each month monthmax(0,*)=[ 0,28,26,30,31,13,31,31,30,31,30,31] ; 2007 monthmax(1,*)=[31,29,30,30,17,30,31,21,30,31,30,31] ; 2008 monthmax(2,*)=[31,28,31,30,31,30,31,31,30,31,30,31] ; 2009 monthmax(3,*)=[31,28,31,30,31,30,31,31,30,26,26,31] ; 2010 monthmax(4,*)=[31,28,31,30,31,18,31,31,30,31,30,31] ; 2011 monthmax(5,*)=[31,29,31,30,31,30,31,31,30,31,30,31] ; 2012 monthmax(6,*)=[31,28,31,30,31,29,31,31,30,31,30,16] ; 2013 monthmax(7,*)=[31,28,31,30,31,30,20,31,30,10,30,31] ; 2014 monthmax(8,*)=[31,28,31,30,31,30,31,31,30,31,30,31] ; 2015 monthmax(9,*)=[31,29,31,30,31,30,31,31,30,31,30,31] ; 2016 monthmax(10,*)=[31,28,31,30,31,30,31,31,30,31,30,31] ; 2017 ; actual number of days in each month monthdays(0,*)= [31,28,31,30,31,30,31,31,30,31,30,31] ; 2007 monthdays(1,*)= [31,29,31,30,31,30,31,31,30,31,30,31] ; 2008 monthdays(2,*)= [31,28,31,30,31,30,31,31,30,31,30,31] ; 2009 monthdays(3,*)= [31,28,31,30,31,30,31,31,30,31,30,31] ; 2010 monthdays(4,*)= [31,28,31,30,31,30,31,31,30,31,30,31] ; 2011 monthdays(5,*)= [31,29,31,30,31,30,31,31,30,31,30,31] ; 2012 monthdays(6,*)= [31,28,31,30,31,30,31,31,30,31,30,31] ; 2013 monthdays(7,*)= [31,28,31,30,31,30,31,31,30,31,30,31] ; 2014 monthdays(8,*)= [ 0, 0, 0, 0, 0,30,31,31,30,31,30,31] ; 2015 monthdays(9,*)= [31,29,31,30,31,30,31,31,30,31,30,31] ; 2016 monthdays(10,*)=[31,28,31,30,31,30, 0, 0, 0, 0, 0, 0] ; 2017 yeardays=total(monthdays,2) imax_all=long(12)*long(12)*long(31)*long(24)*long(60)*long(60) ; x10 = buffer! imax_hr=long(13)*long(12)*long(31)*long(24) imax_d=long(13)*long(12)*long(31) ;*************************************** i_read_SWMF = 1 ; 0 don't read, 1 read i_NoBadMarch = 1 ; 0 means read all data, 1 means skip March 13-15, 2017, when the input was bad i_read_Kyoto = 1 ; 0 = don't read, 1 = read i_doROCstats=1 ; =1 to do the ROC curve statistics, otherwise skip just to make the plot i_ROCnum = 2 ; number of ROC curves to calculate loop_start = 0 ; start of loop for data sets (all available, after-restart gap, and remove bad day) ; probably will be 0 (do all 3) or 2 (do just the "best" data_threshold=[-50.,-30.] ;*************************************** If i_read_SWMF eq 1 then begin ; read in SWMF2011 Dst values Dst_all=fltarr(imax_all) t_Dst_all=dblarr(imax_all) day_Dst_all=intarr(imax_all) Dst_hr_avg=fltarr(imax_hr) Dst_hr_min=fltarr(imax_hr) t_Dst_hr=dblarr(imax_hr) day_Dst_hr=intarr(imax_hr) i_hr=intarr(imax_hr) i_inbadstorm=intarr(imax_hr) Dst_d_avg=fltarr(imax_d) Dst_d_min=fltarr(imax_d) t_Dst_d=dblarr(imax_d) day_Dst_d=intarr(imax_d) i_d=intarr(imax_d) igarb=1 garb='' Dst_read=fltarr(7) iread=long(-1) iread_hr=long(-1) iread_d=long(-1) i_new = 0 run_file='_SWMF2011' run_out='SWMF_2011' i_inbadstorm(*)=0 for iyear=1,3 do begin readyear=iyear+7 ; 2015 and 2016 and 2017 dstyear=2014+iyear ; 2015 and 2016 and 2017 stringyear=STRING(FORMAT='(I4)',dstyear) for imonth=1,12 do begin If monthdays(readyear,imonth-1) gt 0 then begin ; only do months with data If imonth gt 1 then day_offset=total(monthdays(iyear-1,0:imonth-2)) else day_offset = 0 if imonth le 9 then stringmonth='0'+STRING(FORMAT='(I1)',imonth) $ Else stringmonth=STRING(FORMAT='(I2)',imonth) ; for =monthmin(iyear-1,imonth-1),monthmax(iyear-1,imonth-1) do begin ; if iday le 9 then stringday='0'+STRING(FORMAT='(I1)',iday) $ ; Else stringday=STRING(FORMAT='(I2)',iday) openr,1,stringyear+stringmonth+run_file+'_RT_Dst.txt' i_new = 1 For i=1,igarb Do Readf,1,garb while not(EOF(1)) do begin iread=iread+long(1) Readf,1,Dst_read Dst_all(iread)=Dst_read(6) if dst_all(iread) - dst_all(iread) ne 0. then stop if dst_all(iread) lt -1e3 then stop if dst_all(iread) gt 1e3 then stop iday = fix(Dst_read(2)) t_dst_all(iread)= dstyear + ((day_offset + iday -1.) + $ Dst_read(3)/24. + Dst_read(4)/24./60. + Dst_read(5)/24./60./60.)/yeardays(iyear-1) if t_dst_all(iread) - t_dst_all(iread) ne 0. then stop day_dst_all(iread) = iday ; hourly values of SWMF Dst values If i_new eq 1 then begin ; first time through only iread_hr=iread_hr+long(1) this_hr=fix(Dst_read(3)) i_hr(iread_hr)=1 Dst_hr_avg(iread_hr)=Dst_read(6) Dst_hr_min(iread_hr)=Dst_read(6) t_Dst_hr(iread_hr)=dstyear + ((day_offset + iday-1.) + (Dst_read(3)+1.)/24.)/yeardays(iyear-1) day_Dst_hr(iread_hr)=iday endif else if fix(Dst_read(3)) eq this_hr then begin i_hr(iread_hr)=i_hr(iread_hr)+1 Dst_hr_avg(iread_hr)=Dst_hr_avg(iread_hr)+Dst_read(6) If Dst_read(6) lt Dst_hr_min(iread_hr) then Dst_hr_min(iread_hr) = Dst_read(6) endif else begin ; new hr Dst_hr_avg(iread_hr)=Dst_hr_avg(iread_hr)/i_hr(iread_hr) iread_hr=iread_hr+long(1) i_hr(iread_hr) = 1 this_hr = fix(Dst_read(3)) Dst_hr_avg(iread_hr)=Dst_read(6) Dst_hr_min(iread_hr)=Dst_read(6) t_Dst_hr(iread_hr)=dstyear + ((day_offset + iday-1.) + (Dst_read(3)+1.)/24.)/yeardays(iyear-1) day_Dst_hr(iread_hr)=iday endelse ; v3 special if statement: ; Check on whether we are in the bad-data interval of March 13-15, 2017 If i_NoBadMarch eq 1 and dstyear eq 2017 and imonth eq 3 and abs(iday-14) le 1 then begin i_inbadstorm(iread_hr)=1 ; designate as a time in the bad storm interval endif ; Daily values of SWMF Dst values If i_new eq 1 then begin ; first time on a new daily file only i_new = 0 iread_d=iread_d+long(1) this_d=fix(Dst_read(2)) i_d(iread_d)=long(1) Dst_d_avg(iread_d)=Dst_read(6) Dst_d_min(iread_d)=Dst_read(6) t_Dst_d(iread_d)=dstyear + (day_offset + iday-1.)/yeardays(iyear-1) day_Dst_d(iread_d)=iday endif else if fix(Dst_read(2)) eq this_d then begin i_d(iread_d)=i_d(iread_d)+long(1) Dst_d_avg(iread_d)=Dst_d_avg(iread_d)+Dst_read(6) If Dst_read(6) lt Dst_d_min(iread_d) then Dst_d_min(iread_d) = Dst_read(6) endif else begin ; new day Dst_d_avg(iread_d)=Dst_d_avg(iread_d)/i_d(iread_d) iread_d=iread_d+long(1) i_d(iread_d) = 1 this_d = fix(Dst_read(2)) Dst_d_avg(iread_d)=Dst_read(6) Dst_d_min(iread_d)=Dst_read(6) t_Dst_d(iread_d)=dstyear + (day_offset + iday-1.)/yeardays(iyear-1) day_Dst_d(iread_d)=iday endelse endwhile close,1 Dst_d_avg(iread_d)=Dst_d_avg(iread_d)/i_d(iread_d) ; calculate the daily average Dst_hr_avg(iread_hr)=Dst_hr_avg(iread_hr)/i_hr(iread_hr) ; calculate final hourly average of the day ; endfor endif ; check on whether month has data endfor ; month loop endfor ; year loop Dst_all=Dst_all(0:iread) t_Dst_all=t_Dst_all(0:iread) day_Dst_all=day_Dst_all(0:iread) Dst_hr_avg=Dst_hr_avg(0:iread_hr) Dst_hr_min=Dst_hr_min(0:iread_hr) t_Dst_hr=t_Dst_hr(0:iread_hr) day_Dst_hr=day_Dst_hr(0:iread_hr) i_inbadstorm=i_inbadstorm(0:iread_hr) Dst_d_avg=Dst_d_avg(0:iread_d) Dst_d_min=Dst_d_min(0:iread_d) t_Dst_d=t_Dst_d(0:iread_d) day_Dst_d=day_Dst_d(0:iread_d) endif ; end of i_read_SWMF check ; read in Kyoto Dst ;******************************************** If i_read_Kyoto eq 1 then begin Kyoto_dst=fltarr(imax_hr) t_Kyoto_dst=dblarr(imax_hr) day_Kyoto_dst=intarr(imax_hr) Read_Kyoto=dblarr(8) Kyoto_dst_d_avg=fltarr(imax_d) Kyoto_dst_d_min=fltarr(imax_d) t_Kyoto_dst_d=dblarr(imax_d) day_Kyoto_dst_d=intarr(imax_d) i2_d=intarr(imax_d) i_read_2=long(-1) openr,1,'~/Research Codes/SWMF-KyotoDst/DstKyoto.txt' while not(eof(1)) do begin Readf,1,Read_Kyoto if Read_Kyoto(0) ge first_year then begin i_read_2 = i_read_2 + long(1) Kyoto_dst(i_read_2) = Read_Kyoto(7) t_Kyoto_dst(i_read_2) = Read_Kyoto(6) day_Kyoto_dst(i_read_2) = Read_Kyoto(2) ; daily values of Kyoto Dst values If i_read_2 eq 0 then begin ; first time through only iread2_d=long(0) this_d=fix(Read_Kyoto(2)) i2_d(0)=1 Kyoto_Dst_d_avg(iread2_d)=Read_Kyoto(7) Kyoto_Dst_d_min(iread2_d)=Read_Kyoto(7) t_Kyoto_Dst_d(iread2_d)=Read_Kyoto(6) day_Kyoto_Dst_d(iread2_d)=Read_Kyoto(2) endif else if fix(Read_Kyoto(2)) eq this_d then begin i2_d(iread2_d)=i2_d(iread2_d)+1 Kyoto_Dst_d_avg(iread2_d)=Kyoto_Dst_d_avg(iread2_d)+Read_Kyoto(7) If Read_Kyoto(7) lt Kyoto_Dst_d_min(iread2_d) then Kyoto_Dst_d_min(iread2_d) = Read_Kyoto(7) endif else begin ; new day Kyoto_Dst_d_avg(iread2_d)=Kyoto_Dst_d_avg(iread2_d)/i2_d(iread2_d) iread2_d=iread2_d+long(1) i2_d(iread2_d) = 1 this_d=fix(Read_Kyoto(2)) Kyoto_Dst_d_avg(iread2_d)=Read_Kyoto(7) Kyoto_Dst_d_min(iread2_d)=Read_Kyoto(7) t_Kyoto_Dst_d(iread2_d)=Read_Kyoto(6) day_Kyoto_Dst_d(iread2_d)=Read_Kyoto(2) endelse endif endwhile close,1 Kyoto_Dst_d_avg(iread2_d)=Kyoto_Dst_d_avg(iread2_d)/i2_d(iread2_d) ; calculate final daily average Kyoto_dst=Kyoto_dst(0:i_read_2) t_Kyoto_dst=t_Kyoto_dst(0:i_read_2) day_Kyoto_dst=day_Kyoto_dst(0:i_read_2) Kyoto_dst_d_avg=Kyoto_dst_d_avg(0:iread2_d) Kyoto_dst_d_min=Kyoto_dst_d_min(0:iread2_d) t_Kyoto_dst_d=t_Kyoto_dst_d(0:iread2_d) day_Kyoto_dst_d=day_Kyoto_dst_d(0:iread2_d) endif ; i_read_Kyoto check ;******************************************** i_persistence = 0 ; 0 use SWMF, 1 replace SWMF with Kyoto from previous time i_dt_pers = 8 ; hours earlier for the persistence calculations If i_persistence eq 1 then begin ; allow up to 99 hours of persistence shift if i_dt_pers le 9 then stringpers='0'+STRING(FORMAT='(I1)',i_dt_pers) $ Else stringpers=STRING(FORMAT='(I2)',i_dt_pers) run_out = 'Kyoto_pers_'+stringpers+'h' iread=i_read_2 - i_dt_pers ; reset SWMF array length to Kyoto array length, minus persistence shift Dst_all=Kyoto_dst(i_dt_pers:i_read_2) ; shift backwards i_dt_pers time steps t_Dst_all=t_Kyoto_dst(0:iread) ; replace with Kyoto time series day_Dst_all=day_Kyoto_dst(0:iread) ; replace with Kyoto date iread_hr = i_read_2 - i_dt_pers ; hourly values become the shifted Kyoto values Dst_hr_avg=Dst_all Dst_hr_min=Dst_all t_Dst_hr=t_Dst_all day_Dst_hr=day_Dst_all Dst_d_avg=Kyoto_dst_d_avg ; daily values are not used here, just replace with Kyoto values Dst_d_min=Kyoto_dst_d_min t_Dst_d=t_Kyoto_dst_d day_Dst_d=day_Kyoto_dst_d endif ;******************************************** ; search through the two arrays for contingency table yes and no values ; In _ROC version of this code, this is now a loop through the critical value. If i_doROCstats eq 1 then begin delta_t_max = double(1.)/double(365.)/double(24.) ; within one hour in Dst-hr version gap_analysis = 1 ; 0= no check for SWMF restarts, otherwise do the check i_gap = [0,3,3] ; number of hours to search backwards for a gap ;If gap_analysis eq 0 then p_gap = '_yesgap' else p_gap = '_nogap'+STRING(FORMAT='(I1)',i_gap) p_gap='_2gaps' ROC_max = 10 ROC_min = -120 ROC_total = ROC_max - ROC_min If ROC_max le 99 then p_ROC_max = '_ROCmax0'+STRING(FORMAT='(I2)',ROC_max) $ else p_ROC_max = '_ROCmax'+STRING(FORMAT='(I3)',ROC_max) If abs(ROC_min) le 99 then p_ROC_min = '_ROCmin0'+STRING(FORMAT='(I2)',abs(ROC_min)) $ else p_ROC_min = '_ROCmin'+STRING(FORMAT='(I3)',abs(ROC_min)) If i_NoBadMarch eq 1 then p_NoBadMarch = '_NoBadMarch' else p_NoBadMarch = '' fileout=run_out+'_Dst_hr_ROCnSTONE'+p_NoBadMarch+p_gap+p_ROC_max+p_ROC_min TruePos_ROC=fltarr(ROC_total+1,3) FalsePos_ROC=fltarr(ROC_total+1,3) i_good_crit=intarr(ROC_total+1,3) Thresholds=fltarr(ROC_total+1,3) For gap_loop=0,i_ROCnum do begin ; at least twice through, once for STONE and more for ROC calculations print,'Gap loop: ',gap_loop For crit_loop=0,ROC_total do begin ; begin critical value loop critical_value = float(ROC_max-crit_loop) thresholds(crit_loop,gap_loop) = critical_value crit_swmf=critical_value If gap_loop eq 0 then crit_Kyoto=critical_value $ ; this is for STONE calculations else crit_Kyoto=data_threshold(gap_loop-1) ; this is for ROC calculations hits=long(0) misses=long(0) false_a=long(0) corr_neg=long(0) good_Kyoto=fltarr(i_read_2+1) ; here through the next for loop: _hr instead of _d good_swmf=fltarr(i_read_2+1) igood=long(-1) diffsq=0. For i_t1=0,i_read_2 do begin t1=t_Kyoto_dst(i_t1) i_t2=min(where((t_dst_hr + double(.1)*delta_t_max) gt t1)) ; round to nearest time step t2=t_dst_hr(i_t2) If i_t2 ge 0 then begin delta_t = t2 - t1 If gap_analysis eq 0 then begin ; restart check analysis in SWMF results gap_good = 1 endif else begin If i_t2 le i_gap(2) then gap_good = 1 else begin delta_t_swmf = t2 - t_dst_hr(i_t2 - i_gap(2)) If delta_t_swmf gt delta_t_max*(0.5+i_gap(2)) then gap_good = 0 else gap_good = 1 endelse endelse If i_inbadstorm(i_t2) eq 1 then gap_good = 0 ; exclude in the bad storm interval for round 3 if abs(delta_t) le delta_t_max and gap_good eq 1 then begin ; ensure it a usable time igood=igood+long(1) good_Kyoto(igood)=Kyoto_dst(i_t1) good_swmf(igood)=dst_hr_avg(i_t2) If Kyoto_dst(i_t1) le crit_Kyoto then Kyoto_is_crit=1 else Kyoto_is_crit = 0 If dst_hr_avg(i_t2) le crit_swmf then swmf_is_crit=1 else swmf_is_crit = 0 If Kyoto_is_crit eq 1 and swmf_is_crit eq 1 then hits=hits + long(1) If Kyoto_is_crit eq 1 and swmf_is_crit eq 0 then misses = misses + long(1) If Kyoto_is_crit eq 0 and swmf_is_crit eq 1 then false_a = false_a + long(1) If Kyoto_is_crit eq 0 and swmf_is_crit eq 0 then corr_neg = corr_neg + long(1) diffsq=diffsq+(good_Kyoto(igood) - good_swmf(igood))^2 ; If swmf_is_crit eq 1 then print,format='(A25,I10,4F14.5)','SWMF is critical: ', $ ; igood,dst_hr_avg(i_t2),Kyoto_dst(i_t1),t1,t_dst_hr(i_t2) endif endif endfor good_Kyoto=good_Kyoto(0:igood) good_swmf=good_swmf(0:igood) ;******************************************** ; the statistics ;rms_error=sqrt(diffsq/igood) ; igood = N-1 ;good_Kyoto_avg=total(good_Kyoto)/(igood+1) ;good_swmf_avg=total(good_swmf)/(igood+1) ;Kyoto_diffsq=total((good_Kyoto - good_Kyoto_avg)^2) ;swmf_diffsq=total((good_swmf - good_swmf_avg)^2) ;Kyoto_stdev=sqrt(Kyoto_diffsq/igood) ;swmf_stdev=sqrt(swmf_diffsq/igood) ;dst_corr=correlate(good_Kyoto,good_swmf) ;pred_eff = 1.- diffsq/Kyoto_diffsq ;prob_of_det= float(hits)/float(hits+misses) ;prob_of_falsedet = float(false_a)/float(false_a+corr_neg) ;hit_rate = float(hits)/float(hits+misses+false_a) ;hss_numer = 2.*(float(hits*corr_neg) - float(misses*false_a)) ;hss_denom = float(hits+misses)*float(misses+corr_neg) + float(hits+false_a)*float(false_a+corr_neg) ;heidke_skill_score = hss_numer / hss_denom True_Pos = float(hits) / float(hits+misses) False_Pos = float(false_a) / float(false_a+corr_neg) TruePos_ROC(crit_loop,gap_loop)=True_pos FalsePos_ROC(crit_loop,gap_loop)=False_pos i_good_crit(crit_loop,gap_loop)=igood endfor ; end critical value loop print,' Done with that round:',i_good_crit(ROC_total/2,gap_loop),TruePos_ROC(ROC_total/2,gap_loop),FalsePos_ROC(ROC_total/2,gap_loop) endfor ; end gap removal loop endif ; i_doROCstats check ;******************************************** ; Make subsampled arrays for symbol overlay of the lines TruePos_ROC_few=fltarr(ROC_total+1,i_ROCnum+1) FalsePos_ROC_few=fltarr(ROC_total+1,i_ROCnum+1) icount=-1 For i=0,ROC_total,5 do begin ; skip however many you want, or even set it to 1 icount=icount+1 for ROC_loop=0,i_ROCnum do begin TruePos_ROC_few(icount,ROC_loop)=TruePos_ROC(i,ROC_loop) FalsePos_ROC_few(icount,ROC_loop)=FalsePos_ROC(i,ROC_loop) endfor endfor TruePos_ROC_few=TruePos_ROC_few(0:icount,*) FalsePos_ROC_few=FalsePos_ROC_few(0:icount,*) ;******************************************** printfile=1 ; set to 0 if in IDL demo mode If printfile eq 1 then begin ;Print out all of these numbers to a file openw,1,fileout+'.dat' printf,1,'Filename = '+fileout+'.dat' printf,1,'Statistics of '+run_out+' results and real-time Kyoto Dst values: hourly value comparison' printf,1,'Data threshold for ROC calculation (nT): ',data_threshold printf,1,'' printf,1,format='(7A11,A20)','Threshold ',' POD-STONE ', ' POFD-STONE ', $ 'POD-ROC ','POFD-ROC' for crit_loop=0,ROC_total do begin printf,1,format='(F11.2,7F11.6)',float(ROC_max-crit_loop),TruePos_ROC(crit_loop,0),FalsePos_ROC(crit_loop,0), $ TruePos_ROC(crit_loop,1:*),FalsePos_ROC(crit_loop,1:*) endfor close,1 Endif ; printfile check ;************************************************ ; make a plot of the ROC and STONE curves makeplot=1 If makeplot eq 1 then begin !P.multi=[0,0,0] nplots2= 1 ; number of columns xgap=5. xleft=45. ;xgap xright=15. xsize=275. xsz=xleft+xright+nplots2*(xgap+xsize) nplots=1 ; number of rows ygap=5. ysize=150. yupper=25. ylower=40. ysz=ylower+yupper+nplots*(ygap+ysize) y1=ysize/ysz xs1=xgap/xsz xs3=xleft/xsz xs2=xsize/xsz label_x=xs2*.1 ys1=ygap/ysz ys3=ylower/ysz ys2=ysize/ysz label_y=ys2*.05 n1=nplots*nplots2-1 ctk1=.99 ctk2=1. ctk3=.8 !P.charsize=ctk1 Legendlab=['!6STONE curve','ROC Curve (high thr.)','ROC Curve (low thr.)'] ; Plot values For ic=1,2 Do Begin If ic Eq 1 Then Begin set_plot,'X' window,0,title=fileout lth=2 Endif Else Begin set_plot,'ps' device,file=fileout+'.eps',xsize=.01*xsz,ysize=.01*ysz,/inch, $ /port,yoff=1.,scale_factor=1.0,color=8 lth=3 Endelse LOADCT,39 ia=!D.table_size-1 !P.background=ia erase,ia ;lstr=[0,1,2,3,4,5,0] lstr=[0,0,0,0,0,0,0] ;colors=[black,red,blue,orange,green,purple,cyan] lcol=[0,ia-1,ia/4,4*ia/5,3*ia/5,ia/10,2*ia/5] xx=xs3 yy=ys3 !x.range=[0.,1.] !y.range=[0.,1.] plot,FalsePos_ROC(*,0),TruePos_ROC(*,0),xtitle='!6 Probability of False Detection', $ ytitle='!6Prob. of Detection', $ position=[xx,yy,xx+xs2,yy+ys2],charsize=1.,/nodata,thick=lth, $ linestyle=0,title='!6ROC & STONE curves',color=lcol(0),xstyle=1,ystyle=1 oplot,[0.,1.],[0.,1.],linestyle=1,color=lcol(5),thick=lth for gap_loop=0,i_ROCnum do begin oplot,FalsePos_ROC(*,gap_loop),TruePos_ROC(*,gap_loop), $ linestyle=lstr(gap_loop+1),thick=lth,color=lcol(gap_loop+1) oplot,FalsePos_ROC_few(*,gap_loop),TruePos_ROC_few(*,gap_loop),thick=lth,psym=6, $ symsize=.5,color=lcol(1+gap_loop) endfor a=.35 b=.10 c=.025 e=.10 d=.28 f=.025 For i=0,i_ROCnum Do oplot,[a,a+b],[d-i*e,d-i*e],linestyle=lstr(i+1),thick=lth,color=lcol(i+1) for i=0,i_ROCnum Do xyouts,a+b+c,d-i*e-f,Legendlab(i),color=0,charsize=ctk3 If ic NE 1 Then Begin device,/close set_plot,'X' endif endfor endif ; makeplot check ;******************************************** ; binning points for contour overlays valmin=-155. valmax=45. countbinnum = 50 ; sets resolution of the contour plots val_range = valmax - valmin binsize = val_range/float(countbinnum) count_boundaries = findgen(countbinnum+1)*binsize + valmin countcenters = count_boundaries-0.5*binsize countcenters(0)=valmin countcenters(countbinnum)=valmax countbins = intarr(countbinnum+1, countbinnum+1) countbins(*,*) = 0 For i=0,igood do begin If good_kyoto(i) gt valmin then kyoto_bin = min(where(count_boundaries gt good_kyoto(i))) else kyoto_bin=-1 If good_swmf(i) gt valmin then swmf_bin = min(where(count_boundaries gt good_swmf(i))) else swmf_bin=-1 If kyoto_bin*swmf_bin ge 0 then countbins(swmf_bin,kyoto_bin) = countbins(swmf_bin,kyoto_bin)+1 ; If i mod 1000 eq 0 then print,i,good_kyoto(i),kyoto_bin,count_boundaries(kyoto_bin),good_swmf(i),swmf_bin,count_boundaries(swmf_bin),countbins(swmf_bin,kyoto_bin) endfor print,'binned counts: ',max(countbins) binlines = [50,100,150,200,250,300] binlinestyles = [0,0,0,0,0,0] ;[1,2,1,2,1,2] binlinecols = [1,6,4,3,2,5] ;************************************************ ; make a plot of the scatterplot of SWMF and Kyoto For ic=1,2 Do Begin If ic Eq 1 Then Begin set_plot,'X' window,1,title=fileout lth=2 Endif Else Begin set_plot,'ps' device,file=fileout+'_scatter.eps',xsize=.01*xsz,ysize=.01*ysz,/inch, $ /port,yoff=1.,scale_factor=1.0,color=8 lth=3 Endelse LOADCT,39 ia=!D.table_size-1 !P.background=ia erase,ia ;lstr=[0,1,2,3,4,5,0] lstr=[0,0,0,0,0,0,0] ;colors=[black,red,blue,orange,green,purple,cyan] lcol=[0,ia-1,ia/4,4*ia/5,3*ia/5,ia/10,2*ia/5] !x.range=[valmin,valmax] !y.range=[valmin,valmax] plot,good_swmf,good_kyoto,xtitle='!6SWMF Output',ytitle='!6Dst Observations', $ psym=7, symsize=.2,title='!6Dst Values (in nT)',color=0,xstyle=1,ystyle=1 oplot,[valmin,valmax],[valmin,valmax],linestyle=2,color=lcol(5),thick=lth for ROC_loop=0,i_ROCnum-1 do oplot,[-1000.,1000.],[data_threshold(ROC_loop),data_threshold(ROC_loop)], $ linestyle=2,color=lcol(2+ROC_loop),thick=lth contour,countbins,countcenters,countcenters,c_colors=lcol(binlinecols),xstyle=1,ystyle=1, $ levels=binlines,c_linestyle=binlinestyles,/overplot,thick=lth If ic NE 1 Then Begin device,/close set_plot,'X' endif endfor ;************************************************ ; make a plot of the histogram of SWMF and Kyoto beyond some threshold ; binning points for histograms valmin=-155. valmax=45. countbinnum = 20 ; sets resolution of the contour plots val_range = valmax - valmin histbinsize = val_range/float(countbinnum) count_boundaries = findgen(countbinnum+1)*histbinsize + valmin countcenters = count_boundaries-0.5*histbinsize countcenters(0)=valmin countcenters(countbinnum)=valmax skew_vals=fltarr(3) swmf_partial=good_swmf(where(good_kyoto le -30.)) skew_vals(0)=skewness(swmf_partial) swmf_hist30=Histogram(swmf_partial, binsize=histbinsize,max=valmax,min=valmin) swmf_partial=good_swmf(where(good_kyoto le -40.)) skew_vals(1)=skewness(swmf_partial) swmf_hist40=Histogram(swmf_partial, binsize=histbinsize,max=valmax,min=valmin) swmf_partial=good_swmf(where(good_kyoto le -50.)) skew_vals(2)=skewness(swmf_partial) swmf_hist50=Histogram(swmf_partial, binsize=histbinsize,max=valmax,min=valmin) print,'Total,max in SWMF partial range < -30: ',total(swmf_hist30),max(swmf_hist30) p_skew_vals=strarr(3) p_skew_vals(0)='Skew = '+STRING(FORMAT='(F6.2)',skew_vals(0)) p_skew_vals(1)='Skew = '+STRING(FORMAT='(F6.2)',skew_vals(1)) p_skew_vals(2)='Skew = '+STRING(FORMAT='(F6.2)',skew_vals(2)) ctk4=0.7 For ic=1,2 Do Begin If ic Eq 1 Then Begin set_plot,'X' window,2,title=fileout lth=2 Endif Else Begin set_plot,'ps' device,file=fileout+'_histogram_swmf.eps',xsize=.01*xsz,ysize=.01*ysz,/inch, $ /port,yoff=1.,scale_factor=1.0,color=8 lth=3 Endelse LOADCT,39 ia=!D.table_size-1 !P.background=ia erase,ia ;lstr=[0,1,2,3,4,5,0] lstr=[0,0,0,0,0,0,0] ;colors=[black,red,blue,orange,green,purple,cyan] lcol=[0,ia-1,ia/4,4*ia/5,3*ia/5,ia/10,2*ia/5] histcol=[lcol(2),lcol(4),lcol(3)] !x.range=[valmin,valmax] !y.range=[0,410] plot,countcenters,swmf_hist30,xtitle='!6SWMF Output (nT)',ytitle='!6Count', $ title='!6SWMF Histogram',color=0,xstyle=1,ystyle=1,psym=10,thick=lth,/nodata oplot,countcenters,swmf_hist30,color=histcol(2),psym=10,thick=lth oplot,countcenters+0.5,swmf_hist40,color=histcol(1),psym=10,thick=lth oplot,countcenters+1,swmf_hist50,color=histcol(0),psym=10,thick=lth a=valmin+.1*val_range d=.9*410 e=.1*410 for i=0,2 do xyouts,a,d-i*e,p_skew_vals(i),charsize=ctk4,color=histcol(i) If ic NE 1 Then Begin device,/close set_plot,'X' endif endfor skew_vals=fltarr(3) kyoto_partial=good_kyoto(where(good_swmf le -30.)) skew_vals(0)=skewness(kyoto_partial) kyoto_hist30=Histogram(kyoto_partial, binsize=histbinsize,max=valmax,min=valmin) kyoto_partial=good_kyoto(where(good_swmf le -40.)) skew_vals(1)=skewness(kyoto_partial) kyoto_hist40=Histogram(kyoto_partial, binsize=histbinsize,max=valmax,min=valmin) kyoto_partial=good_kyoto(where(good_swmf le -50.)) skew_vals(2)=skewness(kyoto_partial) kyoto_hist50=Histogram(kyoto_partial, binsize=histbinsize,max=valmax,min=valmin) print,'Total,max in Kyoto partial range < -30: ',total(kyoto_hist30),max(kyoto_hist30) p_skew_vals(0)='Skew = '+STRING(FORMAT='(F6.2)',skew_vals(0)) p_skew_vals(1)='Skew = '+STRING(FORMAT='(F6.2)',skew_vals(1)) p_skew_vals(2)='Skew = '+STRING(FORMAT='(F6.2)',skew_vals(2)) For ic=1,2 Do Begin If ic Eq 1 Then Begin set_plot,'X' window,3,title=fileout lth=2 Endif Else Begin set_plot,'ps' device,file=fileout+'_histogram_kyoto.eps',xsize=.01*xsz,ysize=.01*ysz,/inch, $ /port,yoff=1.,scale_factor=1.0,color=8 lth=3 Endelse LOADCT,39 ia=!D.table_size-1 !P.background=ia erase,ia ;lstr=[0,1,2,3,4,5,0] lstr=[0,0,0,0,0,0,0] ;colors=[black,red,blue,orange,green,purple,cyan] lcol=[0,ia-1,ia/4,4*ia/5,3*ia/5,ia/10,2*ia/5] histcol=[lcol(2),lcol(4),lcol(3)] !x.range=[valmin,valmax] !y.range=[0,410] plot,countcenters,kyoto_hist30,xtitle='!6Dst Observations (nT)',ytitle='!6Count', $ title='!6Data Histogram',color=0,xstyle=1,ystyle=1,psym=10,thick=lth,/nodata oplot,countcenters,kyoto_hist30,color=histcol(2),psym=10,thick=lth oplot,countcenters+0.5,kyoto_hist40,color=histcol(1),psym=10,thick=lth oplot,countcenters+1,kyoto_hist50,color=histcol(0),psym=10,thick=lth for i=0,2 do xyouts,a,d-i*e,p_skew_vals(i),charsize=ctk4,color=histcol(i) If ic NE 1 Then Begin device,/close set_plot,'X' endif endfor end