;---------------------------------------------------------------------- ; 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:4,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","winter_2013_14_Stats_1") ; send graphics to PNG file do i_plot = 0,3 data(0,:) = variable_1(0:4,2+i_plot) data(1,:) = variable_2(0:4,2+i_plot) data(2,:) = variable_3(0:4,2+i_plot) data(3,:) = variable_4(0:4,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:4,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.4 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. 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"/) 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 end