; This program makes graphs of TS, POD, FAR, and Bias, as well as a performance plot. It reads four files, but these are in csv format. ; it plots the first three in the performance plot. ; the code could be more elegant, but it works, and it took quite a long time to make it, so I'll stop for now. ; Thanks to Professor Paul Roebber for providing some basic Fortran code that helped me along. Barry Lynn ;---------------------------------------------------------------------- ; This function reads the first line off the CSV file and returns ; the field names. ;---------------------------------------------------------------------- function read_header_from_csv_file(txt_filename) local lines, hdr_split begin lines = asciiread(txt_filename,-1,"string") ;---str_strip removes spaces from begin and end of field names header = str_strip(str_split(lines(0),",")) return(header) end ;---------------------------------------------------------------------- ; This function reads values off a given CSV file, with the header: ; Increment,Day,Time (UTC), Obs, Ens-Mean, E-0, E-1 ,E-2 ,E-3 ,E-4 ,E-5, E-6, E-7 ,E-8 ,E-9 ;---------------------------------------------------------------------- function read_values_from_csv_files(txt_filename,long_name,units) local lines, ntime, nfields, nensem, nt, ne begin delimiter = "," lines = asciiread(txt_filename,-1,"string") print("lines = " + lines) ntime = dimsizes(lines)-1 ; first row is a header header_fields = str_split(lines(1),delimiter) nfields = dimsizes(header_fields) ; nensem = nfields-3 ; first three columns are time ; nensem = nfields-2 ; first two columns are time,threshold nensem = nfields-1 ; first column is time, but this is array value lines_split = str_split_csv(lines(1:),delimiter,0) ; values = tofloat(lines_split(:,3:nensem+2)) ; Start with 4th column of data ; values = tofloat(lines_split(:,2:nensem+1)) ; Start with 3rd column of data values = tofloat(lines_split(:,0:nensem+0)) ; Start with 1stnd column of data return(values) end ; Main code begin plots_collated = new(4,graphic) n_days = 3 ;---------------------------------------------------------------------- ; Read the data off the text files ;---------------------------------------------------------------------- headers_1 = read_header_from_csv_file("prec_stats_day2_12.0km_r12.file") variable_1 = read_values_from_csv_files("prec_stats_day2_12.0km_r12.file","variable_1","in") headers_2 = read_header_from_csv_file("prec_stats_day2_4.0km_r12.file") variable_2 = read_values_from_csv_files("prec_stats_day2_4.0km_r12.file","variable_2","in") headers_3 = read_header_from_csv_file("prec_stats_day2_1.3km_r12.file") variable_3 = read_values_from_csv_files("prec_stats_day2_1.3km_r12.file","variable_3","in") headers_4 = read_header_from_csv_file("prec_stats_day2_12.0km_r12_nocu.file") variable_4 = read_values_from_csv_files("prec_stats_day2_12.0km_r12_nocu.file","variable_1","in") print("headers_1 = " + headers_1) print("headers_2 = " + headers_2) print("headers_3 = " + headers_3) print("variable_1 = " + variable_1) print("variable_2 = " + variable_2) print("variable_3 = " + variable_3) print(dimsizes(variable_1)) ;;;; Make graphs ; xy_2.ncl ; ; Concepts illustrated: ; - Drawing an XY plot with multiple curves ; - Changing the line color for multiple curves in an XY plot ; - Changing the line thickness for multiple curves in an XY plot ; - Drawing XY plot curves with both lines and markers ; - Changing the default markers in an XY plot ; - Making all curves in an XY plot solid ; ;---To plot multiple lines, you must put them into a mulidimensional array. ;data = new((/3,dimsizes(variable_1(:,0))/),float) ;PRINT(DIMSIZES(VARIABLE_1(0:6,0))) data = new((/4,dimsizes(variable_1(0:5,0))/),float) clim = new(4,float) clim(0) = 15.47 clim(1) = 10.91 clim(2) = 6.12 clim(3) = 2.80 wks = gsn_open_wks ("png","Statistics") ; send graphics to PNG file do i_plot = 0,3 data(0,:) = variable_1(0:5,2+i_plot) data(1,:) = variable_2(0:5,2+i_plot) data(2,:) = variable_3(0:5,2+i_plot) data(3,:) = variable_4(0:5,2+i_plot) print("data(1,:) = " + data(0,:)) print("data(2,:) = " + data(1,:)) print("data(3,:) = " + data(2,:)) print("data(4,:) = " + data(3,:)) xaxis = new(dimsizes(data(0,:)-1),float) ;xaxis(0) = 0 xaxis(:) = variable_1(0:5,0) ;---Set plotting parameters res = True ; plot mods desired res@tmXTOn = False res@tmYLFormat = "f" ; remove trailing ".0" ;-------------------------------------------------- ; The time_axis_label function adds additional ; resources to "res" to produce nicely-formatted ; time labels on X axis. This function only works ; if you have a time "units" recognized by the ; cd_calendar function. ;-------------------------------------------------- restick = True res = True ; restick@ttmFormat = "%N/%D %H:%M" ; restick@ttmFormat = "%N/%D" ; time_axis_labels(times_c,res,restick) ;;x_axis_labels(times_c,res,restick) n_t = dimsizes(data(0,:)) ; res@tiMainString = "Disruptive Snowstorms" ; add title res@tiMainString = " " ; set min/max min_h = min(data(:,:)) max_h = max(data(:,:)) ; res@trYMinF = tointeger(min_h) - 50 ; set minimum Y-axis value ; res@trYMaxF = tointeger(max_h) + 50 ; maximum Y-axis value res@trYMinF = 0. res@trYMaxF = 1. if (i_plot.eq.0)then res@trYMinF = 0.2 end if res@trXMinF = 0 ; set minimum X-axis value if (i_plot.eq.3)then res@trYMaxF = 1.5 end if res@trXMinF = 0 ; set minimum X-axis value res@trXMaxF = 2.5 label = (/"12 km"," 4 km","1.3 km"," No-CP"/) colors = (/"blue","red","green","Purple"/) res@xyLabelMode = "Custom" ; label a line res@xyExplicitLabels = label ; text to use res@xyLineLabelFontHeightF = 0.0200 ; font height res@xyLineColors = colors ; label color res@xyLineLabelFontColors = colors ; label color if (i_plot.eq.0)then res@tiYAxisString = "Threat Score" end if if (i_plot.eq.1)then res@tiYAxisString = "Probability of Detection" end if if (i_plot.eq.2)then res@tiYAxisString = "False Alarm Ratio" end if if (i_plot.eq.3)then res@tiYAxisString = "Bias" end if res@tiXAxisString = "Precipitation Threshold (in)" ; ; Similiar resources are xyLineThicknessF and xyLineColor, ; which will affect all lines in the array. ; res@xyLineThicknesses = (/ 3.0, 3.0,3.0,3.0/) res@xyLineColors = (/"blue","red","green","Purple"/) res@tmXBMode = "Explicit" res@tmXBLabelsOn = True res@tmXBValues = xaxis ; location of labels res@tmXBLabels =(/"0.1","0.5","1.0","1.5","2.0","2.5"/) res@gsnDraw = False ; don't draw res@gsnFrame = False ; don't advance frame plot = gsn_csm_xy (wks,xaxis,data,res) ; create plot plots_collated(i_plot) = plot end do resP = True ; panel only resources resP@gsnMaximize = True ; maximize plots gsn_panel(wks,plots_collated,(/2,2/),resP) ; now draw as one plot ; Performance calculations and Graph wks = gsn_open_wks ("png","performance") ; send graphics to PNG file res = True res@xyMarkLineMode = "Markers" ; choose to use markers res@xyMarkers = 16 ; choose type of marker res@trXMinF = 0 ; set minimum X-axis value res@trXMaxF = 1.0 res@trYMinF = 0 ; set minimum X-axis value res@trYMaxF = 1.0 res@tiYAxisString = "POD" res@tiXAxisString = "SR (1-FAR)" dimsx = 6 print(dimsx) xaxis := new((/dimsx*3,3/),float) yaxis := new((/dimsx*3,3/),float) colors := new(dimsx*3,string) markers := new(dimsx*3,integer) mark_sizes := new(dimsx*3,float) ; xaxis(0:dimsizes(SR(0,:) = SR(0,:) ii = 0 do i_var = 0,2 do i = 0,dimsx-1 if (i_var.eq.0)then colors(ii) = "blue" markers(ii) = 3 mark_sizes(ii) = .03 xaxis(0:5,0) = 1.-variable_1(0:5,4) yaxis(0:5,0) = variable_1(0:5,3) end if if (i_var.eq.1)then colors(ii) = "darkgreen" markers(ii) = 12 mark_sizes(ii) = .035 xaxis(6:11,0) = 1.-variable_2(0:5,4) yaxis(6:11,0) = variable_2(0:5,3) end if if (i_var.eq.2)then colors(ii) = "violet" markers(ii) = 15 mark_sizes(ii) = .040 xaxis(12:17,0) = 1.-variable_3(0:5,4) yaxis(12:17,0) = variable_3(0:5,3) end if ii = ii + 1 end do end do ; print(colors) res@xyMarkerColors = colors res@xyMarkers := markers res@xyMarkerSizeF = mark_sizes ; print(res@xyMarkerColors) ; print(xaxis) ; print(yaxis) ; print(x) ; res@tmXBValues := xaxis ; location of labels ; res@tmXBLabels :=(/"0.0","0.1","0.2","0.3","0.4","0.5","0.6","0.7","0.8","0.9","1.0"/) ;print("xaxis = ; yaxis = " + xaxis + " " +yaxis) res@gsnDraw = False ; don't draw res@gsnFrame = False ; don't advance frame plot_stat = gsn_csm_xy (wks,xaxis,yaxis,res) ; create plot ; overlay bias number = 1000 n_pod = number n_sr = number var_ratio = new(1,float) var_double = todouble(var_ratio) bias = new((/n_pod,n_sr/),float) ; csi = new((/n_pod,n_sr/),float) csi = new((/n_pod,n_sr/),integer) dbias = new((/n_pod,n_sr/),float) dcsi = new((/n_pod,n_sr/),integer) ; POD (i_curve) is the variable bias is plotted against (i_point) bias(:,0)=0. a = todouble(1.) b = todouble(1000.) do i_pod=1,n_pod-1 do i_sr = 1,n_sr-1 bias(i_pod,i_sr) = tofloat(i_pod)*.0001/(tofloat(i_sr)*.0001) ; csi(i_pod,i_sr)= 1./(1./(tofloat(i_sr)/1000.)+1./(tofloat(i_pod)/1000.)-1.) var_s = 1.D-3*todouble(i_sr) var_p = 1.D-3*todouble(i_pod) c = a/(a/var_s+a/var_p-a) if (tointeger(c*10000.).le.1001\ .or.tointeger(c*10000).ge.1999.and.tointeger(c*10000).le.2001 \ .or.tointeger(c*10000).ge.2999.and.tointeger(c*10000).le.3001 \ .or.tointeger(c*10000).ge.3999.and.tointeger(c*10000).le.4001 \ .or.tointeger(c*10000).ge.4999.and.tointeger(c*10000).le.5001 \ .or.tointeger(c*10000).ge.5999.and.tointeger(c*10000).le.6001 \ .or.tointeger(c*10000).ge.6999.and.tointeger(c*10000).le.7001 \ .or.tointeger(c*10000).ge.7999.and.tointeger(c*10000).le.8001 \ .or.tointeger(c*10000).ge.8999.and.tointeger(c*10000).le.9001)then ; if (c.le.0.1001) then ;.or.c.ge.19999999.and.c.le.0.2001)then csi(i_pod,i_sr)=tointeger(b*c) else csi(i_pod,i_sr)=csi@_FillValue end if end do end do do i_pod=1,n_pod-1 do i_sr=1,n_sr-1 check = ismissing(csi(i_pod,i_sr)) icsi =0 if (check.eq.False)then if (csi(i_pod,i_sr).eq.100) then icsi=1 end if if (csi(i_pod,i_sr).eq.200) then icsi=2 end if if (csi(i_pod,i_sr).eq.300) then icsi=3 end if if (csi(i_pod,i_sr).eq.400) then icsi=4 end if if (csi(i_pod,i_sr).eq.500) then icsi=5 end if if (csi(i_pod,i_sr).eq.600) then icsi=6 end if if (csi(i_pod,i_sr).eq.700) then icsi=7 end if if (csi(i_pod,i_sr).eq.800) then icsi=8 end if if (csi(i_pod,i_sr).eq.900) then icsi=9 end if end if dcsi(i_pod,i_sr)=icsi ikeep=0 if (bias(i_pod,i_sr).eq.1.0) then ikeep=1 end if if (bias(i_pod,i_sr).eq.1.5) then ikeep=2 end if if (bias(i_pod,i_sr).eq.3.0) then ikeep=3 end if if (bias(i_pod,i_sr).eq.5.0) then ikeep=4 end if if (bias(i_pod,i_sr).eq.10.0) then ikeep=5 end if if (bias(i_pod,i_sr).eq.0.7) then ikeep=6 end if if (bias(i_pod,i_sr).eq.0.5) then ikeep=7 end if if (bias(i_pod,i_sr).eq.0.3) then ikeep=8 end if dbias(i_pod,i_sr)=ikeep end do end do pod_1 = new(n_pod,float) pod_2 = new(n_pod,float) pod_3 = new(n_pod,float) pod_4 = new(n_pod,float) pod_5 = new(n_pod,float) pod_6 = new(n_pod,float) pod_7 = new(n_pod,float) pod_8 = new(n_pod,float) pod_9 = new(n_pod,float) sr_1 = new(n_pod,float) sr_2 = new(n_pod,float) sr_3 = new(n_pod,float) sr_4 = new(n_pod,float) sr_5 = new(n_pod,float) sr_6 = new(n_pod,float) sr_7 = new(n_pod,float) sr_8 = new(n_pod,float) sr_9 = new(n_pod,float) xaxis := new((/10,n_sr/),float) yaxis := new((/10,n_pod/),float) i_val = 0 do i_pod = 1,n_pod-1 do i_sr = 1,n_sr-1 if (dbias(i_pod,i_sr).eq.1)then pod_1(i_val) = tofloat(i_pod)*.001 sr_1(i_val) = tofloat(i_sr)*.001 i_val = i_val + 1 end if end do end do ; an alternative way to find the values, but ii spans the whole range, which is much larger than what I set the possible range ; dcsi1D = ndtooned(dcsi) ;1-d array, same shape as pod_1 and sr_1 ; ii := ind(dcsi1D.eq.1.) ; N := dimsizes(ii) ; newval = fspan(0., (N-1)*0.001, N) ; new array of replacement values ; pod_1(ii) = newval ; sr_1(ii) = newval ; print("dias 2") i_val = 0 do i_pod = 1,n_pod-1 do i_sr = 1,n_sr-1 if (dbias(i_pod,i_sr).eq.2)then pod_2(i_val) = tofloat(i_pod)*.001 sr_2(i_val) = tofloat(i_sr)*.001 i_val = i_val + 1 end if end do end do i_val = 0 do i_pod = 1,n_pod-1 do i_sr = 1,n_sr-1 if (dbias(i_pod,i_sr).eq.3)then pod_3(i_val) = tofloat(i_pod)*.001 sr_3(i_val) = tofloat(i_sr)*.001 i_val = i_val + 1 end if end do end do i_val = 0 do i_pod = 1,n_pod-1 do i_sr = 1,n_sr-1 if (dbias(i_pod,i_sr).eq.4)then pod_4(i_val) = tofloat(i_pod)*.001 sr_4(i_val) = tofloat(i_sr)*.001 i_val = i_val + 1 end if end do end do i_val = 0 do i_pod = 1,n_pod-1 do i_sr = 1,n_sr-1 if (dbias(i_pod,i_sr).eq.5)then pod_5(i_val) = tofloat(i_pod)*.001 sr_5(i_val) = tofloat(i_sr)*.001 i_val = i_val + 1 end if end do end do i_val = 0 do i_pod = 1,n_pod-1 do i_sr = 1,n_sr-1 if (dbias(i_pod,i_sr).eq.6)then pod_6(i_val) = tofloat(i_pod)*.001 sr_6(i_val) = tofloat(i_sr)*.001 i_val = i_val + 1 end if end do end do i_val = 0 do i_pod = 1,n_pod-1 do i_sr = 1,n_sr-1 if (dbias(i_pod,i_sr).eq.7)then pod_7(i_val) = tofloat(i_pod)*.001 sr_7(i_val) = tofloat(i_sr)*.001 i_val = i_val + 1 end if end do end do i_val = 0 do i_pod = 1,n_pod-1 do i_sr = 1,n_sr-1 if (dbias(i_pod,i_sr).eq.8)then pod_8(i_val) = tofloat(i_pod)*.001 sr_8(i_val) = tofloat(i_sr)*.001 i_val = i_val + 1 end if end do end do xaxis(0,:) = sr_1 yaxis(0,:) = pod_1 xaxis(1,:) = sr_2 yaxis(1,:) = pod_2 xaxis(2,:) = sr_3 yaxis(2,:) = pod_3 xaxis(3,:) = sr_4 yaxis(3,:) = pod_4 xaxis(4,:) = sr_5 yaxis(4,:) = pod_5 xaxis(5,:) = sr_6 yaxis(5,:) = pod_6 xaxis(6,:) = sr_7 yaxis(6,:) = pod_7 xaxis(7,:) = sr_8 yaxis(7,:) = pod_8 xaxis(8,:) = sr_9 yaxis(8,:) = pod_9 ; res@xyMarkers = False res@xyDashPattern = 1 label := (/"1","1.5", "3","5","10","0.7","0.5","0.3"/) res@xyLabelMode = "Custom" ; label a line res@xyExplicitLabels := label ; text to use res@xyLineLabelFontHeightF = 0.01250 ; font height res@xyMarkLineMode = "Lines" res@gsnDraw = False ; don't draw res@gsnFrame = False ; don't advance frame plot_bias = gsn_csm_xy (wks,xaxis,yaxis,res) ; create plot pod_1 = new(n_pod,float) pod_2 = new(n_pod,float) pod_3 = new(n_pod,float) pod_4 = new(n_pod,float) pod_5 = new(n_pod,float) pod_6 = new(n_pod,float) pod_7 = new(n_pod,float) pod_8 = new(n_pod,float) pod_9 = new(n_pod,float) sr_1 = new(n_pod,float) sr_2 = new(n_pod,float) sr_3 = new(n_pod,float) sr_4 = new(n_pod,float) sr_5 = new(n_pod,float) sr_6 = new(n_pod,float) sr_7 = new(n_pod,float) sr_8 = new(n_pod,float) sr_9 = new(n_pod,float) xaxis := new((/10,n_sr/),float) yaxis := new((/10,n_pod/),float) i_val = 0 do i_pod = 1,n_pod-1,3 do i_sr = 1,n_sr-1,3 if (i_val.le.n_pod-1)then if (dcsi(i_pod,i_sr).eq.1)then pod_1(i_val) = tofloat(i_pod)*.001 sr_1(i_val) = tofloat(i_sr)*.001 i_val = i_val + 1 end if end if end do end do i_val = 0 do i_pod = 1,n_pod-1,3 do i_sr = 1,n_sr-1,3 if (i_val.le.n_pod-1)then if (dcsi(i_pod,i_sr).eq.2)then pod_2(i_val) = tofloat(i_pod)*.001 sr_2(i_val) = tofloat(i_sr)*.001 i_val = i_val + 1 end if end if end do end do i_val = 0 do i_pod = 1,n_pod-1,3 do i_sr = 1,n_sr-1,3 if (i_val.le.n_pod-1)then if (dcsi(i_pod,i_sr).eq.3)then pod_3(i_val) = tofloat(i_pod)*.001 sr_3(i_val) = tofloat(i_sr)*.001 i_val = i_val + 1 end if end if end do end do i_val = 0 do i_pod = 1,n_pod-1,2 do i_sr = 1,n_sr-1,2 if (i_val.le.n_pod-1)then if (dcsi(i_pod,i_sr).eq.4)then pod_4(i_val) = tofloat(i_pod)*.001 sr_4(i_val) = tofloat(i_sr)*.001 i_val = i_val + 1 end if end if end do end do i_val = 0 do i_pod = 1,n_pod-1,2 do i_sr = 1,n_sr-1,2 if (i_val.le.n_pod-1)then if (dcsi(i_pod,i_sr).eq.5)then pod_5(i_val) = tofloat(i_pod)*.001 sr_5(i_val) = tofloat(i_sr)*.001 i_val = i_val + 1 end if end if end do end do i_val = 0 do i_pod = 1,n_pod-1,2 do i_sr = 1,n_sr-1,2 if (i_val.le.n_pod-1)then if (dcsi(i_pod,i_sr).eq.6)then pod_6(i_val) = tofloat(i_pod)*.001 sr_6(i_val) = tofloat(i_sr)*.001 i_val = i_val + 1 end if end if end do end do i_val = 0 do i_pod = 1,n_pod-1 do i_sr = 1,n_sr-1 if (i_val.le.n_pod-1)then if (dcsi(i_pod,i_sr).eq.7)then pod_7(i_val) = tofloat(i_pod)*.001 sr_7(i_val) = tofloat(i_sr)*.001 i_val = i_val + 1 end if end if end do end do i_val = 0 do i_pod = 1,n_pod-1 do i_sr = 1,n_sr-1 if (i_val.le.n_pod-1)then if (dcsi(i_pod,i_sr).eq.8)then pod_8(i_val) = tofloat(i_pod)*.001 sr_8(i_val) = tofloat(i_sr)*.001 i_val = i_val + 1 end if end if end do end do i_val = 0 do i_pod = 1,n_pod-1 do i_sr = 1,n_sr-1 if (i_val.le.n_pod-1)then if (dcsi(i_pod,i_sr).eq.9)then pod_9(i_val) = tofloat(i_pod)*.001 sr_9(i_val) = tofloat(i_sr)*.001 i_val = i_val + 1 end if end if end do end do xaxis(0,:) = sr_1 yaxis(0,:) = pod_1 xaxis(1,:) = sr_2 yaxis(1,:) = pod_2 xaxis(2,:) = sr_3 yaxis(2,:) = pod_3 xaxis(3,:) = sr_4 yaxis(3,:) = pod_4 xaxis(4,:) = sr_5 yaxis(4,:) = pod_5 xaxis(5,:) = sr_6 yaxis(5,:) = pod_6 xaxis(6,:) = sr_7 yaxis(6,:) = pod_7 xaxis(7,:) = sr_8 yaxis(7,:) = pod_8 xaxis(8,:) = sr_9 yaxis(8,:) = pod_9 ; res@xyMarkers := 1 res@xyDashPattern = 0 label := (/"0.1","0.2", "0.3","0.4","0.5","0.6","0.7","0.8","0.9"/) res@xyLabelMode = "Custom" ; label a line res@xyExplicitLabels := label ; text to use res@xyLineLabelFontHeightF = 0.01250 ; font height ; res@xyMarkers := markers res@gsnDraw = False ; don't draw res@gsnFrame = False ; don't advance frame colors := (/"black","black","black","Purple"/) plot_csi = gsn_csm_xy (wks,xaxis,yaxis,res) ; create plot res@gsnDraw = True ; don't draw res@gsnFrame = True ; don't advance frame overlay(plot_stat,plot_csi) overlay(plot_stat,plot_bias) draw(plot_stat) frame(wks) ; print(x) end